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Executive  Summary 


Classification  and  Discrimination  of  Sources  with  Time-Varying 
Frequency  and  Spatial  Spectra 


ONR,  Grant  no.  N000 14-98- 1-0176 
Moeness  Amin  (PI) 

This  report  presents  the  results  of  the  research  performed  under  ONR  grant 
number  N00014-98-1-0176  over  the  period  of  October  1st,  02  to  September  30th,  03.  The 
research  team  working  on  this  project  consists  of:  Prof.  Moeness  Amin  (PI),  Prof.  Yimin 
Zhang  (Research  Professor),  Dr.  Gordon  Frazer  (Visiting  Research  Professor  from 
DSTO,  Australia),  Mr.  Weifeng  Mu,  Mr.  Baha  Obeidat,  and  Mr.  Behzad  Mohammadi 
(Graduate  Students).  We  have  also  collaborated  with  Prof.  X-G.  Xia  from  University  of 
Delaware,  Prof.  H.  Ge  from  NJIT,  and  Prof.  A.  Zoubir  from  Darmstadt  University  of 
Technology,  Germany. 

The  research  objectives  in  FY03  evolved  around  the  detection  of  complex 
Doppler  target  signatures  using  time-frequency  signal  representations.  In  additions  to 
achieving  these  objectives,  the  research  has  progressed  on  two  other  fronts;  namely, 
Polarimetric  Sensor  Array  Processing  and  Subband  Array  Implementations  for  Space- 
Time  Processing.  Four  book  chapters,  five  journal  article,  and  ten  conference  papers  have 
been  generated  from  our  research  efforts  in  FY03.  Below  is  the  summary  of  the  research 
accomplishments  in  each  of  the  above  three  areas.  Copies  of  the  principle  publications 
are  included  in  the  report. 

1.  Over  the  Horizon  Radar 

We  have  developed  a  high-resolution  time-varying  signature  estimation  of  weak  target 
signals  for  Over-the-Horizon  Radars  (OTHRs).  A  significant  problem  in  OTHR  is  robust  high- 
resolution  Doppler  processing  of  maneuvering  targets.  The  complex  Doppler  signatures  present  in 
these  cases  smear  in  Doppler  and  have  reduced  detectability  and  localization.  We  have  introduced 
a  powerful  clutter  suppression  technique  combined  with  high-resolution  time-varying  signature 
estimation  method,  specifically,  in  low  signal-to-clutter  ratio  environment.  This  technique  is 
based  on  an  integration  of  autoregressive  filtering,  adaptive  chirp/chirplet  transform,  and  adaptive 
bilinear  transform.  The  new  technique  is  most  applicable  to  scenarios  where  resolution  of 
multipath  returns  may  clearly  reveal  accelerating/decelerating  characteristics  of  maneuvering 
targets  which  prove  significant  for  target  classifications.  In  these  situations,  the  Fourier  transform 
and  many  existing  time-frequency  methods  often  fail. 

2.  Polarimetric  Sensor  Array  Processing 

We  have  developed  polarimetric  time-frequency  distributions  for  sensor  arrays.  A  new 
technique  that  encompasses  the  target  time-varying  Doppler  signature  as  well  as  the  target 
polarization  signature  for  applications  to  high-resolution  imaging  and  direction  finding  has  been 
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introduced.  Both  signatures  are  estimated  and  processed  by  multi-antenna  receivers.  The  signal 
polarization  information  empowers  the  recently  introduced  spatial  time-frequency  distributions, 
leading  to  improved  source  discrimination  and  angle-of-arrival  estimation.  The  proposed 
technique  is  a  general  approach  that  utilizes  diversity  and  distinctions  amongst  targets  based  on 
their  polarizations,  directions,  and  time-frequency  signal  representations.  Coupling  of  those  three 
source  characteristics  has  proven  superior  to  dealing  with  each  one  separately. 

3.  Subband  Array  Implementations  for  Space-Time  Processing 

Space-time  adaptive  processing  (STAP)  systems  have  been  shown  to  be  effective  in 
suppressing  undesired  signals  in  Radar  and  mobile  communication  applications.  The  high 
complexity  and  slow  convergence,  however,  are  often  the  hurdles  in  practical  implementations  of 
STAP  systems.  Subband  array  implementations  provide  near-optimal  steady  state  performance 
with  reduced  implementation  complexity  and  improved  convergence  performance.  We  continued 
to  analyze  the  performance  of  several  subband  array  processing  schemes.  In  particular,  we  have 
evaluated  the  effect  of  decimation  on  the  mean  square  errors  and  examined  the  advantages  of 
performing  subband  processing  on  the  received  input  data  as  well  as  the  desired  (reference) 
signal. 
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High  Resolution  Time-Frequency 
Distributions  for  Maneuvering  Target 
Detection  in  Over-the-Horizon  Radars 


Yimin  Zhang  and  Moeness  G.  Amin 

Center  for  Advanced  Communications,  Villanova  University,  Villanova,  PA  19085,  USA 

Gordon  J.  Frazer 

ISR  Division,  DSTO,  Edinburgh,  Australia 


Abstract 

A  novel  high-resolution  time-frequency  representation  method  is  proposed  for  source  detection  and 
classification  in  over-the-horizon  radar  (OTHR)  systems.  A  data-dependent  kernel  is  applied  in  the 
ambiguity  domain  to  capture  the  target  signal  components,  which  are  then  resolved  using  the  root-MUSIC 
based  coherent  spectrum  estimation.  This  two-step  procedure  is  particularly  effective  for  analyzing  a 
multi-component  signal  with  time- varying  complex  time-Doppler  signatures.  By  using  the  different  time- 
Doppler  signatures,  important  target  maneuvering  information,  which  is  difficult  to  extract  using  other 
linear  and  bilinear  time-frequency  representation  methods,  can  be  easily  revealed  using  the  proposed 
method. 

This  work  was  supported  by  the  Office  of  Naval  Research  under  Grant  NG0014-98-1-0176. 
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I.  Introduction 


By  exploiting  the  reflective  and  refractive  nature  of  high-frequency  (HF)  radiowave 
propagation  through  the  ionosphere,  over-the-horizon  radars  (OTHRs)  perform  wide-area 
surveillance  at  long  range  well  beyond  the  limit  of  the  horizon  of  conventional  line-of-sight 
(LOS)  radars.  OTHR  systems  have  become  an  important  tool  for  wide-area  surveillance 
[1],  [2],  [3],  [4],  [5], 

A  significant  problem  in  OTHR  is  robust  high-resolution  Doppler  processing  of  accel¬ 
erating  or  decelerating  targets.  This  arises  during  aircraft  and  ship  target  maneuver  and 
during  observations  of  rockets  during  boost  phase  and  mid-course  flight.  The  complex 
Doppler  signatures  present  in  these  cases  reveal  important  information  about  the  target. 

Most  OTHR  systems  use  classical  Doppler  processing  where  one  Doppler  spectrum 
is  computed  using  one  full  coherent  integration  time  (CIT,  typically  1  -  100  seconds 
in  OTHR).  Some  systems  use  overlapped  Doppler  processing  to  provide  a  spectrogram 
analysis  of  time- varying  Doppler.  Accelerating/decelerating  targets  smear  in  Doppler  and 
have  reduced  detectability  and  localization.  The  smearing  reduces  resolution  and  can 
obscure  important  multi-component  Doppler  features. 

There  are  numerous  time-frequency  distributions  (TFDs)  other  than  the  spectrogram 
[6],  [7],  [8].  Many  TFDs  provide  superior  localization  in  time  and  Doppler  frequency. 
Previous  applications  of  time-frequency  signal  representations  to  OTHR,  however,  have 
generally  been  disappointing.  The  fundamental  challenge  with  OTHR  is  that  the  TFD 
must  retain  its  desirable  resolution  and  concentration  properties  in  the  presence  of  clutter 
that  is  typically  40dB  or  more  stronger  than  the  target  (although  possibly  localized  in  a 
different  region  of  time- Doppler). 

The  objective  of  this  paper  is  to  investigate  and  extend  recent  developments  in  data- 
dependent  TFDs  to  the  problem  of  robust  high-resolution  analysis  of  time-varying  OTHR 
target  returns.  Of  particular  interest  is  the  problem  of  multi-component  target  signal  de¬ 
tection  and  identification  where  important  information  concerning  the  maneuvering  target 
can  be  revealed.  Such  information  is  significant  for  the  classification  of  maneuvering  tar¬ 
gets. 

This  paper  is  organized  as  follows.  In  Section  II,  the  signal  model  is  introduced  and 
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the  Doppler  characteristics  of  the  target  signal  are  investigated.  Section  III  provides 
the  clutter  suppression  method  used  in  this  paper.  In  Section  IV,  a  high-resolution  time- 
Doppler  processing  method  is  proposed  and  applied  to  the  underlying  OTHR  applications. 
Section  V  shows  some  more  computational  results. 

II.  Signal  Model 

A.  Signal  Model 

Fig.  1(a)  illustrates  an  stylised  OTHR  system.  For  simplicity  of  mathematical  analysis, 
we  adopt  the  flat  ground  model,  as  shown  in  Fig.  1(b). 

The  received  signal,  after  pulse  or  sweep  matched  filtering  and  beamforming  at  the 
receiver,  is  expressed  as 

y(t)  =  x(t)  +  u(t)  +  n(t),  (1) 

where  u(t)  is  the  clutter,  and  x(t)  is  the  return  signal  from  the  target,  expressed  as 

x(t)  =  Ae~j^dt+dr)/c,  (2) 

where  A  is  a  complex  scalar  representing  the  propagation  loss  and  phase,  dt  and  dr  are  the 
respective  one-way  slant  range  between  the  transmitter  and  the  target  and  between  the 
target  and  the  receiver,  c  denotes  the  speed  of  light,  and  uc  =  2i r/c  is  the  carrier  radian 
frequency.  In  (1),  n(t )  is  internal  and  external  noise,  whose  power  is  small  in  a  typical 
situation  with  strong  target  signal-to-noise  ratio.  Therefore,  the  noise  term  is  ignored  in 
this  paper. 

In  a  typical  OTHR  scenario,  as  shown  in  Fig.  1(a),  in  addition  to  the  path  directly 
reflected  from  the  ionosphere,  there  is  reflection  at  the  ground  near  the  target.  Denote 
li  and  1%  as  the  propagation  distance  of  the  two  paths,  respectively,  and  dt  and  dr  as  the 
respective  one-way  slant  range  between  the  transmitter  and  the  target  and  between  the 
target  and  the  receiver,  respectively.  Then,  dt  takes  value  of  either  k  and  l2  and  so  does 
dr}  The  corresponding  path  losses  will  be  denoted  as  A\  and  A2.  Therefore,  the  received 

^ome  OTHR  systems  are  bistatic  in  which  case  the  transmitter  and  receiver  are  in  different  locations.  In  this 
case,  the  range  from  the  target  to  the  transmitter  and  that  from  the  target  to  the  receiver  are  different.  We  have 
ignored  this  difference,  however,  as  it  does  not  significantly  affect  our  results. 
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signal  consists  of  four  combination  paths  which  result  in  the  following  three  multipath 
components: 

x(t)  =  Aie-jWc2h/c  +  A2e~]UJc2h/c  +  A3e-julc{h+h)/c.  (3) 

We  will  refer  to  the  path  (/i  :  h)  as  path  I,  path  (/2  :  h)  as  path  II,  and  the  combination 
of  path  (/i :  l2)  and  (/2  :  h)  as  path  III. 

Based  on  the  flat  ground  model  approximation  illustrated  in  Fig.  1(b),  the  slant  ranges 
l\  and  I2,  respectively,  can  be  expressed  in  terms  of  the  range  distance  R,  the  ionosphere 
height  H,  and  the  aircraft  height  h,  as 

h  =  (R2  +  (2H  -  h)2)l/\  l2  =  {R2  + (2H  +  h)2)l/2 .  (4) 


To  clearly  reveal  the  relationship  between  the  slant  ranges  and  the  parameters,  we  take 
into  account  the  fact  that  R  >>  H  >■  h  hold  for  a  typical  scenario.  Then,  the  above 


expressions  can  be  approximated  as 


2  H2  +  2  Hh 


B.  Doppler  Characteristics 


The  flight  of  an  aircraft,  in  general,  consists  of  horizontal  and  elevation  movements.  In 
this  section,  we  consider  the  Doppler  frequency  characteristics  of  the  aircraft’s  movement 
in  the  two  different  dimensions. 


As  an  aircraft  flies,  R ,  and  possibly  also  h,  become  functions  of  t.  The  height  of 
ionosphere,  H,  is  also  slowly  time-varying.  However,  we  assume  that  H  is  a  constant  over 
the  processing  time  interval.  From  (5),  we  obtain 

^ d2f~mvM+ (6) 

where  K{t)  —  (1  —  2 H2/R2(t)),  vR(t)  —  dR(t)/dt  is  the  target  velocity  in  the  range 
direction  toward  the  radar,  and  vc(t)  —  dh(t)/dt  is  the  ascending  velocity  of  the  target. 


The  Doppler  frequencies  of  the  three  different  paths  are  then  obtained  as 

2fcdk(t)  2/c  4  fcH 

m = -~dT K  -  mfcVc{)’ 


fn(t)  = 


2/c  dl2(t) 


-K(t)vR(t) 


4 fJ1_ 

R(t)c 


Vc(t), 


fm{t)  = 


fcdh(t)  +  dl2(t) 


:-K(t)vR(t). 
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From  the  above  discussion,  it  is  evident  that,  while  the  dominant  Doppler  component 
2 fcK(t)vR(t)/c  is  shared  by  all  the  three  paths  and  reveals  the  target  velocity  in  the 
range  direction,  the  small  Doppler  difference  between  the  paths  is  a  function  of  h(t)  and, 
therefore,  reveals  another  important  information  on  how  the  target  moves  in  the  elevation 
direction. 

In  this  paper,  we  consider  an  often  encountered  scenario  of  a  maneuvering  aircraft  as 
an  example.  In  this  case,  the  target  makes  a  180°  turn  in  T  =  30.72  seconds  to  change 
height  and  direction.  This  time  interval  corresponds  to  6  revisits  (blocks),  and  each  block 
contains  256  samples  (slow  time  samples  from  the  radar).  The  parameters  used  in  the 
analysis  and  simulations  are  listed  in  Table  1.  All  the  multipath  signals  are  considered  to 
fall  within  the  same  range  cell. 

The  range  is  expressed  as 

R(t)  =  R(Q)  -  VR^T  sin  (|0  (8) 

and  the  height  is  expressed  as 

h(t)  =  h( 0)  +  Vc®T  1  -  cos  J  .  (9) 

The  cross-range  movement  is  not  considered  because  it  does  not  significantly  contribute  to 
Doppler  frequency.  Substituting  (8)  and  (9)  to  (5)  yields  the  following  Doppler  frequency 

f(t)  =  fR(t)+  fh(t),  (10) 


where  /r  (f)  is  the  Doppler  frequency  caused  by  the  change  of  R(t)  and  is  expressed  as 

/«(»)  =  cos(^V),  (11) 


and  fh(t)  in  (10)  is  the  Doppler  frequency  caused  by  the  change  of  height  h(t)  and  is 
expressed,  for  path  I,  as 


4 vc(t)fcH  . 


It  is  easy  to  confirm  that  fh,n{t )  =  — )  and  fh,ni{t)  —  0  for  all  t. 

The  time-Doppler  signatures  are  plotted  in  Fig.  2.  The  dominant  Doppler  component 
is  proportional  to  the  target  velocity  in  the  slant  range  direction,  and  the  small  Doppler 
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difference  between  the  three  paths  is  proportional  to  the  ascending  velocity  of  the  target. 
This  difference  provides  important  information  on  how  the  target  moves  in  the  eleva¬ 
tion  direction.  The  maximum  one-side  Doppler  difference  corresponding  to  the  maximum 
ascending  speed  1500  m/min  =  25  m/s  is  about  1.17  Hz. 

The  frequency  resolution  in  the  underlying  system  is  A /  =  fs/N  =  50/256  =  0.195  Hz. 
The  detection  of  such  a  small  Doppler  difference  is  possible  through  the  application  of  dis¬ 
crete  Fourier  transform  (DFT)  to  the  received  signal,  provided  that  v(t)  and,  subsequently, 
2 fcv(t)/c,  are  fixed  over  the  CIT  of  256  samples.  However,  when  vr  is  not  a  constant, 
which  occurs  if  the  target  is  accelerating  or  decelerating,  ascending  or  descending,  or, 
changing  its  direction,  2 fcv(t)/c  becomes  time- varying  and  the  conventional  DFT-based 
approach  does  not  provide  high  Doppler  resolution  even  with  a  long  CIT  [10].  In  this  case, 
the  detection  and  estimation  of  the  Doppler  shift  caused  by  a  change  of  h  become  difficult. 
The  presence  of  strong  clutter  adds  more  difficulties  to  the  underlying  problem. 

Conventional  methods  based  on  spectrogram  and  other  TFDs  smear  the  target’s  Doppler 
signature  and  cannot  provide  satisfactory  resolution  performance.  The  smearing  reduces 
resolution  and  is  likely  to  obscure  important  multi-component  time-Doppler  signatures.  To 
realize  high-resolution  Doppler  detection  and  estimation,  we  must  first  proceed  with  clutter 
suppression  followed  by  an  effective  time-Doppler  processing  method.  The  combination  of 
the  two  methods  clearly  reveals  the  interested  time-Doppler  signatures. 

III.  Clutter  Suppression 

We  consider  TFD  methods  to  achieve  high-resolution  time-Doppler  localization.  In  the 
underlying  problem,  TFDs  are  referred  to  as  time-Doppler  distributions  (TDDs).  The 
most  commonly  used  TFD  is  the  Wigner-Ville  Distribution  (WVD).  The  WVD  of  signal 
y(t)  is  defined  as 

Wyy{t,  /)  =  J  y(t  +  r/2)y*(t  -  r/2)e"Ja’Tdr  (13) 

where  the  superscript  denotes  complex  conjugate.  All  integrals  without  limits  imply 
integration  from  — oo  to  +00.  Substituting  (1)  in  (13),  the  WVD  of  y(t)  can  be  written 
in  terms  of 

Wyy(t,  f)  =  Wxx(t ,  /)  +  Wuu(t,  f)  +  Wxu(t,  f )  +  Wux(t,  /),  (14) 


6 


where  the  first  two  terms  are,  respectively,  the  autoterms  of  the  target  signal  and  the 
clutter,  and  the  other  two  are  their  crossterms. 

In  a  typical  OTHR  receiver,  the  clutter  is  much  stronger  (typically  30  to  50  dB)  than 
the  target  signal.  Without  substantial  suppression  of  the  clutter,  the  TDD  autoterm  of 
the  target  will  be  significantly  obscured  by  the  clutter  autoterm  as  well  as  the  crossterms 
between  the  clutter  and  signal. 

Clutter  often  has  high  correlation  to  that  at  its  neighboring  range  cells  and  cross-range 
cells.  Based  on  this  property,  clutter  mitigation  methods  using  received  signals  from  other 
range  and  cross-range  cells  have  been  proposed  in,  for  example,  [9],  [10].  However,  in  this 
paper,  the  received  signal  from  other  range  and  cross-range  cells  are  not  used. 

We  point  to  the  fact  that  the  clutter  is  highly  localized  in  low  frequencies  and  can  be 
well  modeled  as  an  autoregressive  (AR)  process  [11],  [12].  Therefore,  the  clutter  can  be 
substantially  suppressed  by  using  the  AR  pre- whitening  techniques.  Denote  P  as  the  order 
of  the  AR  model,  the  AR  polynomial  parameters  a(t),t  =  0,  ...,P  can  be  estimated  via 
the  modified  covariance  method  [13]. 

Filtering  the  received  signal  y(t )  through  a  finite  impulse  filter  (FIR),  constructed  using 
the  AR  polynomial  parameters,  results  in  the  pre-whitened  signal: 

z(t)  =  y(t)  *  a(t )  =  x(t )  *  a(t )  +  u(t)  *  a(t)  A  zx(t)  +  zu(t),  (15) 

where  “*”  denotes  the  convolution  operator. 

In  this  paper,  the  target  signal  calculated  in  Section  II  is  overlaid  to  real  OTHR  clutter 
data.  We  assume  that  Ai  —  A2,  and  A$  =  Ai  +  A2  =  2 A\.  The  order  of  the  AR  model 
should  be  chosen  to  maximize  the  signal- to-clutter  ratio  (SCR).  The  order  of  the  AR  model 
is  set  to  a  unit  value  (P— 1).  The  spectrogram  of  block  3,  which  corresponds  to  the  256 
samples  from  10.24  to  15.36  seconds,  before  and  after  the  AR  pre-whitening  is  shown  in 
Fig.  3.  It  is  seen  that,  while  the  clutter  is  substantially  suppressed  by  more  than  40  dB, 
the  target  signal  is  only  partially  affected  when  its  Doppler  frequency  is  very  close  to  that 
of  the  clutter.  Figs.  4  and  5  show  the  WVDs  of  the  y(t)  and  z(t)  before  and  after  the 
pre-whitening.  The  WVDs  are  computed  from  the  interpolated  data  sequence  to  show  the 
full  Doppler-frequency  range.  While  it  is  often  difficult  to  identify  the  target  in  the  WVD 
before  pre- whitening  (Fig.  4),  the  target  signature  can  now  be  somewhat  identified  in  Fig. 
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5.  Further  and  key  improvement  in  resolutions  of  the  target  signature  components  can  be 
achieved  by  using  the  techniques  highlighted  in  Section  4. 

IV.  High-Resolution  Time-Doppler  Processing 

Even  after  substantial  clutter  suppression,  the  result  in  Fig.  5  does  not  reveal  clear  time- 
Doppler  signatures.  There  are  numerous  TFDs  other  than  the  spectrogram  and  the  WVD 
which  provide  superior  localization.  Previous  applications  of  time-frequency  signal  rep¬ 
resentation  techniques  to  the  OTHR  problem,  however,  has  generally  been  disappointing 
because  problem  is  particularly  difficult  and  demanding. 

To  achieve  chirp  signal  detection,  discrimination,  and  classification,  we  propose  time- 
Doppler  estimation  based  on  adaptive  kernel  and  high-resolution  time-Doppler  localiza¬ 
tion.  Bilinear  TDDs  as  well  as  their  Fourier  transforms  (i.e,  ambiguity  functions  and  local 
auto-correlation  functions)  are  considered. 

A.  Time-Doppler  Distributions  and  Adaptive  Kernel 

In  the  following,  we  assume  that  each  component  of  the  return  signal  from  the  target 
can  be  approximated  as  a  chirp  over  the  period  of  one  block,  i.e., 

x(t)  =  J2Aiejiait+Pxt2/2)  (16) 

1  =  1 

is  considered  for  a  time  period.  Such  approximation  permits  us  to  obtain  important  signal 
information,  as  discussed  in  Section  II,  from  the  received  data  signal. 

To  estimate  the  chirp  rate  of  the  signal,  it  is  common  to  examine  the  ambiguity  function. 
The  ambiguity  function  of  z(t)  is  defined  as 

Az(0,  r)  —  J  z(t  +  rj2)z*(t  —  rj2)e^tdt  (17) 

t 

where  6  and  r  are,  respectively,  the  frequency-lag  and  time-lag  variables.  Similar  to  the 
TDD,  the  ambiguity  function  can  be  decomposed  into  two  autoterms  and  two  crossterms. 
One  important  property  of  the  ambiguity  function  is  that  all  signal  autoterms  pass  through 
the  origin,  whereas  the  crossterms  are  often  away  from  the  origin.  For  a  multi-component 
parallel  chirp  signal,  the  ambiguity  function  shows  linear  signatures  depending  on  the  sig¬ 
nal  chirp  rate.  Therefore,  unlike  the  time-Doppler  domain,  in  which  a  two-dimensional 
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search  is  required,  the  chirp  rate  in  the  ambiguity  domain  can  be  estimated  by  a  one¬ 
dimensional  search.  The  reducution  in  computations  make  the  ambiguity-domain  attrac¬ 
tive  for  chirp  rate  estimation. 

The  chirp  rate  can  be  estimated  by  searching  for  the  peaks  of  the  following  Q  function 
[14] 

<5(0  =  J \Az(r  cos£,rsin£)|dr.  (18) 

In  the  case  considered,  peaks  possibly  appear  at  =  —1/  tan_1(/?x)  and  £„  =  —1/  tan_1(/?u), 
where  j3x  and  fiu  are  the  chirp  rates  of  the  signal  and  the  principle  component  of  the  resid¬ 
ual  clutter,  respectively. 

Based  on  the  chirp  rate  estimation,  an  adaptive  kernel  can  be  designed.  We  construct 
a  kernel  whose  passband  only  captures  the  target  signal  chirp  rate.  The  clutter  will  be, 
subsequently,  mitigated  in  the  ambiguity  domain  due  to  its  distinct  orientation  compared 
to  the  target  signal.  For  an  estimated  chirp  rate  £x,  the  following  adaptive  kernel  is 
constructed  to  encompass  the  autoterm  ambiguity  function  of  the  target  signal,  i.e., 

<t>a{e,T)  =  e-<mu/-2  (19) 

where  a  is  the  kernel  width,  and 

d2(0,r)  =  92  +  r2  —  (0sin|x  +  tcos£x)2.  (20) 

The  adaptive  kernel  suppresses  the  clutter  and  noise,  as  well  as  all  crossterms. 

The  adaptive  chirp  TDD  is 

c,(t, »)4EE t)M0, (21) 

/7r  e  t 

The  above  distribution  has  substantially  suppressed  clutter  and  noise,  as  well  as  the 
crossterm  TDDs  between  the  multi-component  signal.  The  adaptive  TDD  is  shown  in 
Fig.  7  for  the  received  signal  at  block  3. 

B.  High  Resolution  Time-Doppler  Localization 

In  [14],  chirp  MUSIC  was  introduced  for  the  estimation  of  the  Doppler  frequencies  at 
each  time  index  t.  The  estimated  auto-correlation  function  Rx(t,  r)  is  used  to  construct  a 
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data  matrix  for  MUSIC  spectrum  estimation.  However,  the  resulting  matrix  is,  in  general, 
not  positive  definite.  Therefore,  in  [14],  the  filtered  ambiguity  function  is  transformed 
to  the  time-frequency  domain,  and  only  the  positive  part  of  the  TFD  is  considered  for 
the  construction  of  the  auto-correlation  function.  This  method,  although  showing  good 
time-Doppler  localization  in  high  signal-to-noise  ratio  (SNR)  situations,  is  computationally 
inefficient  because  spectrum  estimation  is  implemented  for  each  time  index.  In  addition, 
the  estimated  time-Doppler  signature  is  not  always  consistent  with  the  true  values,  par¬ 
ticularly  in  low  SNR  scenarios.  Therefore,  it  is  not  a  candidate  for  application  in  the 
underlying  OTHR  applications. 

In  this  paper,  we  obtain  the  auto-correlation  directly  from  the  filtered  ambiguity  function 
as 

Rx(t,r)  =  -^  J  Ax(0,T)<f>a(O,T)e~jetdd.  (22) 

Because  signal  components  with  single  chirp  rate  are  involved,  the  auto-corelation  function 
Rx,i(t ,  t)  of  each  chirp  component  has  the  form 

Rx^r)  =  A2ie^ai+^\  (23) 

which  is  dependent  of  t.  Such  dependence  can  be  removed  by  using  the  estimated  value, 
fix  =  —1/  tan(£r).  From  Rx,i(t,  r),  the  time-independent  auto-correlation  function  is  esti¬ 
mated  as 

Rx(t)  =  J  Rx(t,r)e~j^xtTdt.  (24) 

The  coherent  integration  yields  coherent  MUSIC  subspace  estimation  of  at’s  for  improved 
performance.  The  vector  Rx(r)  is  considered  as  raw  data  sequence,  rather  than  as  co- 
variance  elements  adopted  in  [14],  to  ensure  the  positive  definiteness  of  the  covariance 
matrix  for  spectrum  estimation.  In  our  simulations,  the  root-MUSIC  algorithm  is  used  for 
computational  convenience.  Only  one  root-MUSIC  operation  is  required  for  each  block. 
The  chirp  signatures  at  different  times  are  then  constructed  using  the  estimated  chirp  rate 
and  ccj’s. 

In  Fig.  8,  the  coherent  time- varying  root-MUSIC  spectrum  is  shown  for  block  3.  Despite 
the  low  SCR,  the  time-Doppler  signatures,  along  with  the  Doppler  frequency  difference 
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information,  are  estimated  clearly  and  consistently.  Simulation  results  for  all  other  blocks 
also  confirmed  successful  Doppler  signature  estimation. 

C.  TDD  Magnitude  Compression  and  More  Simulation  Results 

The  existence  of  strong  time- Doppler,  value  at  some  discrete  points,  however,  may  some¬ 
time  create  undesired  time-Doppler  signatures.  Because  the  desired  TDDs  typically  show 
much  more  consistent  signature  over  all  samples  with  its  true  chirp  rate,  we  propose  the 
use  of  the  following  magnitude  compression  of  Cx(t,u)), 

C'x(t,c a)  =  \Cx(t,  w)|7  sign[Cx(t,  o>)],  (25) 

where  0  <  7  <  1.  Our  experience  suggests  that  7  should  take  value  between  0.1  and  0.5. 
Cx(t,  10)  is  used  instead  to  estimate  the  auto-correlation  function  Rx(t,r )  in  (26). 

When  the  TDD  magnitude  compression  is  performed,  the  local  auto-correlation  function 
Rx(t,r )  should  be  obtained  from 

R,(t,r)=±-  ECi((,»y.  (26) 

U) 

In  Fig.  9,  the  time-Doppler  signatures  obtained  from  the  proposed  method  is  shown  for 
the  entire  turning  process  of  the  aircraft.  In  the  computations,  the  results  are  calculated 
from  six  blocks,  each  of  256  samples.  7  =  0.2  are  used  for  each  block.  The  theoretical 
values  of  the  Doppler  signatures  are  overlaid  in  the  plot.  It  is  evident  that  the  proposed 
method  provides  stable  and  consistent  estimation  of  the  Doppler  signatures  over  different 
situations. 

To  show  the  importance  of  applying  magnitude  compression,  we  plotted  in  Fig.  10  the 
time-Doppler  signatures  obtained  without  the  magnitude  compression.  It  is  seen  that, 
while  most  time-Doppler  signatures  are  correctly  estimated,  one  component  in  block  4  is 
not.  The  reason  is  simply  that,  in  the  process  of  clutter  suppression,  signal  component 
with  close  spectrum  to  the  clutter  may  loss  part  of  its  signal  power. 

V.  Conclusions 

In  this  paper,  a  novel  method  has  been  proposed  for  high-resolution  time-Doppler  sig¬ 
nature  localization  applied  to  over-the-horizon  radar  systems.  By  combining  AR  pre¬ 
whitening  for  effective  clutter  suppression,  time-frequency  based  signal  discrimination, 
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and  coherent  high-resolution  spectrum  analysis,  the  proposed  method  provides  robust  es¬ 
timation  of  time- varying  Doppler  signature  in  low  signal-to-clutter  ratio  (SCR)  scenarios. 
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TABLE  I 


Major  Parameters 


Parameter 

Notation 

Value 

initial  range 

m 

2000  km 

ionosphere  height 

H 

350  km 

aircraft  initial  height 

Ml o) 

10000  m 

maximum  range  speed 

VR 

500  km/hr 

maximum  climbing  speed 

Vc 

1500  m/min 

carrier  frequency 

fc 

20  MHz 

waveform  repetition  frequency 

fs 

50  Hz 

samples  per  block 

N 

256  samples 
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Fig.  1.  (a)  OTHR  systems,  (b)  Flat  ground  approximation. 
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Fig.  6.  Q  function  calculated  (block  3). 
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Fig.  7.  Adaptive  time-Doppler  distribution  of  the  pre- whitened  received  signal  (block  3) 
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Fig.  8.  Estimated  time-Doppler  signature  via  chirp  root-MUSIC  algorithm  (block  3). 


Fig.  9.  Estimated  time-Doppler  signature  of  all  blocks  with  magnitude  compression  (‘+’  marks  show  the 
theoretical  Doppler  frequencies  for  the  three  paths) . 
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Fig.  10.  Estimated  time-Doppler  signature  of  all  blocks  without  magnitude  compression  (‘-h’  marks  show 
the  theoretical  Doppler  frequencies  for  the  three  paths). 
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Abstract 

In  over-the-horizon  radar  (OTHR)  systems,  the  signal -to-clutter  ratio  (SCR)  used  for  moving 
target  detection  is  very  low.  For  slowly  moving  targets  such  as  ships,  the  SCR  is  typically  from  - 
50  dB  to  -60  dB  and  their  Doppler  frequencies  are  close  to  that  of  the  clutter.  For  maneuvering 
targets,  such  as  aircrafts  and  missiles,  the  Doppler  frequencies  are  time-varying  when  a  long 
integration  time  is  considered.  When  a  target  does  not  move  uniformly,  the  Fourier  transform 
based  target  detection  techniques,  including  super-resolution  spectrum  techniques,  may  fail  to 
work  appropriately.  In  such  situations,  the  Doppler  signatures  are  time-varying  and,  therefore, 
time-frequency  analysis  techniques  can  be  used  for  maneuvering  target  detection.  In  addition, 
clutter  rejection  is  also  required  for  target  detection  due  to  the  low  SCR.  The  existing  adaptive 
clutter  rejection  algorithms  combine  clutter  rejection  with  spectrum  analysis  methods  which 
usually  assume  uniformly  moving  target  (i.e.,  sinusoidal  Doppler  signature)  models.  In  this  paper, 
we  propose  an  adaptive  clutter  reject  algorithm  together  with  the  adaptive  chirplet  transform 
technique  for  maneuvering  target  detection  in  a  multipath  environment.  Simulation  results  using 
simulated  maneuvering  target  signal  with  received  raw  OTHR  clutter  data  show  that  targets  with 
SCR  below  -50  dB  can  be  detected  by  using  the  proposed  algorithm. 
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1.  Introduction 


Over-the-horizon  radar  (OTHR)  systems  have  been  widely  used  to  detect  and  track  targets,  such  as 
aircrafts  and  surface  ships,  in  wide-area  surveillance  at  long  ranges  [1-5].  The  existing  OTHR  detection 
and  tracking  algorithms  are  based  on  the  assumption  that  the  Doppler  frequency  of  each  target  is  constant 
or  at  least  approximately  constant  during  each  dwell.  Targets  are  detected  from  amplitude  peaks  away 
from  the  zero  frequency.  The  detection  capability  of  an  algorithm  depends  on  the  SCR  and  the  Doppler 
resolution  which,  in  turn,  depends  on  the  length  of  the  coherent  integration  time  (CIT).  In  the  existing 
Fourier  transform  based  techniques  for  maneuvering  target  detection  and  tracking,  there  is  a  trade-off 
between  the  CIT  length,  SCR,  and  the  Doppler  resolution.  For  a  slowly  or  uniformly  moving  target,  such 
as  a  ship,  the  Fourier  transform  based  techniques  work  well,  where  a  long  CIT  can  be  used  for  clutter 
spread  reduction.  However,  for  a  fast  maneuvering  target,  such  as  fast  boat,  aircraft,  or  missile,  in  Fourier 
transform  based  techniques,  a  long  CIT  can  not  be  used  and,  therefore,  the  Doppler  resolution  is  limited. 
In  such  situations,  time-frequency  analysis  becomes  an  important  technique  for  effective  maneuvering 
target  detection  and  tracking.  Time-frequency  analysis  methods  [6-9]  have  found  wide  applications  in 
radar,  [10-15].  Because  the  radar  return  signals  from  maneuvering  targets  have  chirp-like  characteristics, 
a  new  Doppler  processing  method  based  on  the  adaptive  chirplet  transform  (ACT),  instead  the  Fourier 
transform,  is  proposed  in  this  paper.  With  the  adaptive  chirplet  transform  technique,  the  CIT  can  be 
substantially  extended  and,  therefore,  the  Doppler  resolution  can  be  improved  compared  with  the  Fourier 
transform  based  techniques. 

In  an  OTHR  system,  the  detection  of  slow  targets  is  often  difficult,  due  to  clutter  from  the  ground 
or  ocean  and  the  low  SCR  (typically  about  -50  dB  to  -60  dB).  Therefore,  clutter  rejection  is  necessary  to 
improve  the  target  detect  capability.  The  available  clutter  rejection  algorithms  include  Fourier  transform 
based  adaptive  clutter  rejection  method  recently  proposed  by  Root  [2],  and  super-resolution  spectrum 
estimation  algorithms,  for  example,  [16,17],  In  this  paper,  an  adaptive  clutter  rejection  algorithm  is 
proposed.  After  clutter  rejection,  the  ACT  is  then  applied  to  the  clutter-mitigated  signal,  which  makes  the 
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energy  of  the  maneuvering  targets  concentrated.  By  using  the  proposed  algorithm,  the  maneuvering 
targets  can  be  correctly  detected  even  when  the  SCR  is  below  -50  dB. 

This  paper  is  organized  as  follows.  In  Section  2,  the  OTHR  signal  model  is  briefly  reviewed.  In 
Section  3,  the  adaptive  chirplet  transform  for  OTHR  is  introduced.  In  Section  4,  the  adaptive  clutter 
rejection  algorithm  is  proposed.  In  Section  5,  simulation  results  of  maneuvering  target  detection  are 
provided. 


2.  OTHR  Signal  Model  and  Problem  Description 

In  this  section,  we  first  describe  the  OTHR  signal  model  and  the  conventional  OTHR  processing  for 
uniformly  moving  targets,  and  then  the  problem  of  interest  in  this  paper  for  maneuvering  targets. 

2.1  OTHR  Signal  Model  for  OTHR  Processing 

After  the  low  pass  filtering  and  sampling,  the  received  signal  s(m,n )  for  a  target  p  with  ground 
range  r  is  (see  for  example  [4,  5]) 

+  CD 

where  n ,  m  ,T0,  (Of ,  dp,  f  ,  B ,  Tc  and  £nm  are  the  fast  time  sample  index,  chirp  pulse  index,  the 

minimum  delay,  Doppler  frequency  shift,  the  two-way  slant  rang,  waveform  repetition  frequency, 
bandwidth  of  radar,  coherent  integration  time  and  additive  noise,  respectively.  Ap  in  equation  (1)  is  the 

amplitude  of  the  received  signal  from  target  (source)  or  clutter  p .  The  power  Ap  of  the  signal  has  the 

P  G  (7 

expression  as  [21  ]A2  = — - — --x - —  x  Ae ,  where  R,  and  Rr  are  the  distances  (meters)  of  transmitter 

'  AnRt2  4xRr2  ‘ 

and  receiver  of  radar  to  the  target  respectively,  Pt  is  the  power  (watts)  of  antenna  with  gain  Gt ,  o  is  the 
cross  section  (in  square  meters)  of  the  target,  and  Ae  is  the  effective  area  of  antenna  aperture.  From  (1), 
we  find  that  the  signal  part  in  s(m,n)  in  terms  of  index  n  is  a  complex  sinusoidal  signal.  It  is  also  a 


s(m, n)  =  Ap  exp (jcopmTc ) cxp|^ 


dP  1 

<op-2nBf{-?--T0)\nTs 

c  J 
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sinusoidal  function  of  index  m  if  the  Doppler  frequency  cop  does  not  change  with  m  .  In  this  case,  a  two- 
dimensional  discrete  Fourier  transform  over  m  and  n  provides  the  range-Doppler  surface  S(m',ri) .  The 
received  clutter  signal  is  the  signal  of  coming  from  the  ground  and  surface  of  sea.  Signal  of  clutter  is 
spread  in  frequency,  it  is  not  just  appear  in  zero  frequency.  In  the  real  OTH  radar,  there  are  two  high  peak 
correspond  to  the  Bragg  lines  away  from  zero  frequency  which  make  the  target  detection  more  difficult. 
For  a  particular  OTHR  processing  algorithm,  the  target  detection  capability  depends  on  the  SCR. 
Therefore,  in  order  to  improve  the  target  detection  performance,  one  can  increase  the  range,  Doppler 

c 

resolution,  and/or  the  SCR.  The  range  resolution,  A r  -  —  ,  depends  on  the  bandwidth  B  of  radar. 

2  B 

2k 

However,  the  Doppler  resolution,  A 0)  -  — ,  depends  on  the  CIT  Tc ,  which  is  chosen  at  the  receiver. 

Targets  and  clutter  with  Doppler  difference  less  than  Aco  are  located  in  one  Doppler  cell.  One  Doppler 
cell  is  divided  into  k  smaller  cells  and  the  SCR  is  then  increased  by  k  times  if  the  CIT  increases  k  times. 
The  assumption  here  is  that  the  target  moves  uniformly  within  the  CIT  interval.  This  assumption, 
however,  may  not  hold  when  the  CIT  is  long. 


2.2  Problem  Description  on  Maneuvering  Target  Detection 

For  a  maneuvering  target,  the  signal  Doppler  frequency  a>p  in  (1),  due  to  the  target  motion  or  non- 

uniform  motion  of  electron  density  distribution  in  ionospheric  [17],  is  no  longer  constant  but  time 
varying.  Consider  a  moving  target  with  initial  velocity  v  and  acceleration  a  in  the  direction  of  slant 

4k  4k 

range.  The  Doppler  frequency  cop in  (1)  is  cop{t)  -—(v  +  at)  .  The  Doppler  spread  is  Aa>p  =  —aTc, 

A  A 

A  (0  2 aT2 

and  thus,  the  number  of  Doppler  cells  that  the  target  energy  spreads  over  is - —  = - — .  Therefore, 

A  CD  X 

when  the  target  moves  uniformly,  i.e.,  a  -  0 ,  the  target  energy  is  always  concentrated  in  a  single  Doppler 


cell.  It  becomes  different,  however,  when  the  target  moves  non-uniformly,  i.e.,  a  0 .  As  an  example,  let 
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us  assume  2a!  X  =  \ .  In  this  case,  the  target  energy  spreads  over  Tc2  Doppler  cells.  This  implies  that,  if 
the  CIT  Tc  increases  k  times,  the  number  of  Doppler  cells  over  which  the  target  energy  spreads 

increases  k 2  times.  Therefore,  in  this  case,  the  SCR  in  Doppler  reduces  k 2  times  compared  to  that  in  the 
uniform  moving  target  case.  This  tells  us  that,  for  a  maneuvering  target,  the  CIT  increase  does  not  benefit 
the  OTHR  target  detection  if  the  Fourier  transform  is  used  in  the  Doppler  processing.  We  next  want  to 
propose  an  adaptive  chirplet  transform  (ACT)  in  the  Doppler  processing  that  may  take  the  advantage  of 
the  long  CIT  no  matter  whether  the  target  moves  uniformly  or  not. 

3.  Chirp  Signal  Detection  and  Adaptive  Chirplet  Transform 

In  OTHR,  the  received  signal  corresponding  to  a  range  cell  is  typically  a  multi-component  signal 
with  time-varying  frequency  signatures  corresponding  to  the  multiple  targets  with  different  velocities  and 
the  clutter.  In  this  section,  we  first  review  the  Wigner-Ville  distribution  (WVD)  and  Radon-Wigner 
transform  (RWT)  for  multiple  chirp  detection.  We  then  describe  the  adaptive  chirplet  transform  method 
which  is  used  in  the  simulation  in  Section  5. 

3.1  Chirp  Signal  Detection  Using  Radon-Wigner  Transform 

WVD  of  a  signal  s(t)  is  defined  as  follows  [7,8] 

Ws{t,0))  =  J  s(t  +  ^)s\t-^)exp(-j(OT)dT,  (2) 

—oo  ^  " 

where  variables  t  and  0)  represent  the  time  and  frequency,  respectively.  WVD  has  the  highest  resolution 
for  a  single  chirp  signal,  but  its  major  disadvantage  is  the  presence  of  artificial  cross-terms  caused  by  the 
quadratic  multiplication  nature.  For  a  signal  containing  multiple  linear  chirps,  the  desired  WVD  auto¬ 
terms  are  straight  lines  in  the  Wigner  plane,  while  the  undesired  cross-terms  are  manifested  as  the  high 
frequency  oscillating  characteristics.  To  suppress  the  cross-terms,  we  consider  the  Radon-Wigner 
transform  (RWT)  which  takes  advantage  of  the  above  oscillating  properties  by  integrating  the  WVD 
along  lines  with  different  chirp  rate  and  frequency  shift  combinations.  A  large  part  of  the  WVD  cross- 
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terms  is  cancelled  to  each  other  through  the  integration,  and  the  residual  part  of  the  cross-terms  can  be 
further  reduced  in  the  Radon-Wigner  plane  by  noting  the  fact  that  the  RWT  auto-  and  cross-terms  have 
different  characteristics.  Therefore,  a  proper  mask  can  be  applied  to  the  RWT  to  reduce  the  cross-terms 
with  minimum  distortion  to  the  auto-terms.  The  WVD  with  substantially  suppressed  cross-terms  can  be 
obtained  by  transforming  the  masked  RWT  back  to  the  Wigner  plane  [20].  It  is  proved  in  [20]  that  the 
WVD  auto-terms  after  the  inverse  Radon  transform  of  the  masked  coefficients  are  the  same  as  those  in 
the  original  WVD.  These  WVD  auto-terms  are  then  used  to  estimate  the  instantaneous  Doppler 
frequencies  of  targets.  For  other  instantaneous  Doppler  frequency  estimation  methods,  see  for  example 
[6]. 

For  multi-component  signals  with  approximately  equal  magnitudes,  the  RWT  filtering  in  the 
Radon-Wigner  plane  is  effective.  However,  when  the  magnitudes  of  the  signal  components  are 
significantly  different,  the  method  may  not  be  effective  because  the  cross-terms  between  stronger  signal 
components  may  be  larger  than  the  auto-terms  of  weaker  components.  In  this  case,  a  weaker  signal  may 
be  shaded  in  the  presence  of  strong  cross-terms  and  can  hardly  be  detected.  In  OTHR  systems,  the  signal 
echo  from  small  targets,  such  as  small  boats,  are  often  much  weaker  than  that  of  the  clutter  even  after  the 
clutter  cancellation.  In  this  case,  the  method  we  introduced  above  can  be  used  to  detect  the  strongest 
signal  component  and  then  remove  it  from  the  original  signal.  This  procedure  is  repeated  until  all  the 
signal  components  are  detected. 

3. 2  Adaptive  Chirplet  Transform  for  High  Order  Time-Varying  Frequency  Signals 

For  a  long  CIT,  the  received  signal  from  a  maneuvering  target  is  no  longer  a  linear  chirp  signal. 
When  the  time-varying  frequency  is  a  higher  order  polynomial  of  time,  the  signal  can  be  expressed  as  a 
combination  of  several  linear  chirps  over  different  time  intervals.  Such  kind  of  signal  representations, 
introduced  by  Mann  and  Haykin  [15],  is  called  chirplet  transform. 
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The  procedure  of  chirplet  decomposition  of  a  signal  is  first  to  estimate  the  chirp  rates  a, ,  a2, 


aNo  of  s(t)  over  different  segments,  and  the  respective  chirps  w, (r)  =  exp jj(— ocf  )j  are  then 

constructed.  For  a  given  frame  {hk(t),keZ}  and  N0  chirp  rates,  a  new  chirplet  frame  \fik(t)Uj(t)Jc,ie  Z} 
is  obtained.  Based  on  this  chirplet  frame  {hk (r)«,(r)},  s(t)  is  divided  as 

5(f)  =  ZZca^A(0K,(0  (3) 

1=1  k 

where  Cik  =  (s(t),hk(Lt)ui(t))  are  the  frame  decomposition  and  fyik(t),kez}  is  the  dual  frame  of 
{hk  (t),  ke  Z},  (•,•)  represents  the  inner  product,  g,  are  arbitrary  weights  satisfying 

£ft=l.  0<?;<1.  (4) 

i= 1 

For  details  about  (3),  see  [13].  To  have  an  efficient  frame  decomposition,  {hk( t),  keZj  should  include 


functions  with  different  time  and  frequency  bandwidths  and  center  (mean)  locations.  For  example,  the 


following  modulated  Gaussian  functions 


M 0  =  f— 1  exp [-nO-hf  +  J  A+^r(t~tk) | 

K  J  l  2 


are  usually  used,  where  yk ,  <pk  are  parameters  that  control  the  envelop  and  phase  of  the  chirplet,  and  fik 
and  tk  denote  the  frequency  and  time  centers,  respectively. 

Next,  we  consider  how  to  construct  the  chirplet  frame.  Radon-Wigner  transform  can  be  used  to 
estimate  these  chirp  rates.  For  a  given  signal  s(t) ,  chirp  rate  ax  is  obtained  by  searching  the  largest  peak 
in  the  Radon-Wigner  plane  after  taking  the  RWT  of  the  signal.  We  then  obtain  the  chirplet  frame 

(X 

{hk (r)u,(/)}  by  modulating  the  frame  {hk(t)}  in  (5)  with  «1(f)  =  exp(y'-^-r2).  We  next  estimate  which 
element  in  the  modified  frame  {hk  (t)ul  (r)}  optimally  matches  the  signal  s(t)  and  denote  this  element  as 


«j  (t)hk  (t)  where 
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K  (0  =  arg  min  il  s(t)  -  {t)K  (t)  ||  l 

k  [  llnt(t)ll  J 


Define  signal  ^  (t)  as 

(s(t),ux(t)hk(t)) 

5,  (r)  =  s(t) - - - n,  ( t)hk  (t). 

1  ll/z,(t)ll  1  *'  (?) 

By  repeating  the  same  procedure  of  s(t)  to  s,(t) ,  we  obtain  the  chirp  rate  a2  corresponding  to  the 


second  largest  component  of  s(t) .  Let  u2(t)  =  exp (j-^-t2) ,  we  obtain 


.  ..  .  f n  .  .  (s,(t),u2(t)hk(t))u2(t)hk(t) 

*2  *  \  1  ihkm 


.  ,  <S,(t),M2(0\(0>  ...  .. 

s2(t)  =  s,(t) - - - u2(t)hk  (t). 

2  1  11/1,(011  2  *2  <9) 

Repeating  the  above  procedure,  all  signal  components  can  be  obtained,  and  signal  s(t)  can  be  expressed 

as  s(t)  =  Sj :(/).  Based  on  the  above  decomposition,  the  instantaneous  frequencies  of  all  signal 

I 

components  can  be  obtained  and  then  used  for  OTHR  target  detection. 

One  can  see  that  the  search  in  (6)  and  (8)  is  four-dimensional  and  thus  has  a  high  computational 
complexity.  In  order  to  reduce  the  complexity,  a  sub-optimal  adaptive  chirplet  transform  algorithm  with 
two  dimensional  search  is  given  in  [13]  and  is  summarized  as  follows: 

Step  1.  Chirp  rate  aj ,  and  frequency  6\  are  estimated  by 

(or, ,  <y, )  =  arg  max^  (a,  co)\ ,  (10) 

(a, to) 

where  Ds(a,a> )  is  the  RWT  of  s{t) . 

Step  2.  Let  g(t)  denote  a  match  filter  with  center  frequency  0)x .  Also,  denote  ux(t)  =  exp(y^-r2)and 


hi(t)  =  g(t)°(s(t)u1(t)) 
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where  °  is  the  convolution  operator. 

Step  3.  The  coefficient  C, ,  in  (3)  is  obtained  as 

C, ,  =(s(t),h1(t)ul(t))=  js(t)hl(t)u*(t)dt  . 


Step  4.  Let 


sl(t)=Cl  lhl(t)  Ui(t),  and  y,(t)  =  s(/)-.s,(f)  . 


Step  5.  Set  s  =  y,  (?) . 

Step  6.  Stop  if  the  energy  of  s(t)  is  small  enough,  otherwise  go  to  Step  1. 
Other  adaptive  chirplet  transforms  can  be  found  in,  for  example  [9]. 


(12) 


(13) 


4.  Clutter  Rejection 

In  OTHR,  clutter  is  a  multi-component  signal  with  much  stronger  power  than  that  of  target  signal. 
To  achieve  effective  target  detection,  clutter  rejection  is  necessary  before  ACT  is  applied.  A  clutter 
rejection  algorithm  using  adaptive  Fourier  transform  was  proposed  by  Root  [2].  There  are  some  other 
algorithms  [17,18],  in  which  adaptive  clutter  rejection  and  maximum  likelihood  target  detection  are 
combined  based  on  the  sinusoidal  target  signal  model.  In  this  section,  adaptive  clutter  rejection  algorithms 
are  discussed,  which  are  independent  of  target  detection  methods  and  do  not  assume  any  target  signal 
model.  After  clutter  rejection,  time-frequency  analysis  can  be  used  to  make  the  energy  of  maneuvering 
target  focused. 

4.1  “Adaptive  Noise  Canceling”  Method  Used  to  Clutter  Rejection 

To  effectively  suppress  the  clutter,  we  notice  the  fact  that  the  clutter  has  high  space  correlation  to 
its  neighboring  range  cells.  The  correlation  coefficient  may  be  as  high  as  0.8-0.9  [17,  18].  Therefore,  the 
received  signals  at  neighboring  range  cells  can  be  used  to  estimate  the  clutter  covariance  matrix  of  the 
current  range  cell.  The  idea  of  adaptive  noise  canceller  can  be  used  to  the  underlying  OTHR  clutter 
rejection  problem.  An  adaptive  noise  canceller  is  a  dual-input,  closed-loop  adaptive  feedback  system  [19], 
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which  makes  use  of  the  signal  d(n)  received  at  the  primary  sensor  and  the  signal  v,(n)  received  at  a 
reference  sensor.  The  signal  received  d(n)  at  the  primary  sensor  is  composed  of  the  interested  signal 
s(n)and  additive  noise  v0(n) ,  i.e., 

d(n)  =  s(n)  +  v0(n )  (14) 

It  is  assumed  that  the  signal  s(n)  and  noise  v0(n)  are  uncorrelated  to  each  other,  and  v,(n)  is  correlated 
to  the  noise  v0(n)  but  is  uncorrelated  to  the  signal  s(n) .  The  reference  signal  v,(n)  is  used  to  estimated 
the  noise  component  in  d(n) , 

y(n)  =  ''jTwk(n)vl(n-k)  (15) 

*=o 

where  wk(n)  are  the  adjustable  tap  weights  of  the  adaptive  filter.  The  filter  output  y(n)  is  subtracted 
from  the  primary  signal  d(n) ,  resulting  in  the  following  error  signal 

e(n )  =  d(n)  -  y(n)  =  s(n)  +  v0(«)  -  y(n) .  (16) 

The  error  signal  is  used  to  adjust  the  tap  weights  of  the  adaptive  filter.  The  error  signal  e(n)  is  the  overall 
system  output,  which  contains  the  desired  signal  s(n)  with  the  noise  suppressed. 

In  the  OTHR  target  detection,  the  desired  signal  is  the  echo  from  targets,  whereas  the  undesired 
signal  is  clutter  and  noise.  As  the  OTHR  signal  is  usually  processed  after  beamforming,  a  reference 
sensor  receiving  only  the  clutter  is  not  available.  However,  as  we  mentioned  before,  the  clutter  in  one 
range  cell  has  high  correlation  with  that  of  its  neighboring  range  cells  whereas  the  target  signals  do  not. 
Therefore,  the  signals  received  at  neighboring  range  cells  can  be  used  for  clutter  suppression,  resemble  to 
the  signal  received  at  the  reference  sensor  in  an  adaptive  noise  canceller. 

4.2  Adaptive  Clutter  Rejection  Algorithm 

In  this  subsection,  an  adaptive  signal  subspace  method  is  used  to  clutter  rejection.  Let  sc(t)  be  the 
received  signal  after  range  compression  in  the  current  interested  range  cell,  sni(t),...,snfi(t)  are  received 
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signals  after  range  compression  from  N  neighboring  range  cells.  Define  the  following  signal  vectors 
Sc  and  5,  constructed  for  the  current  and  the  ith  neighboring  range  cells  from  the  all  M  samples  over  the 


Se=[sA0),s  -l)]T . 


Then,  the  covariance  matrix  of  clutter  and  external  noise  can  be  estimated  by 


mrSisr 
r=—n — 
mr 


where 


is  the  correlation  coefficient  between  the  received  signal  vectors  at  the  current  range  cell  and  the  ith 
neighboring  cell,  and  y  is  a  positive  scalar,  typically  takes  value  between  1  and  2. 

The  introduction  of  term  |^,|r  in  (19)  allows  weighting  differently  the  contributions  of  the 

neighboring  range  cells  depending  on  their  respective  correlation  coefficients  with  the  current  range  cell. 
By  doing  this,  the  contributions  from  less  correlated  range  cells  can  be  effectively  eliminated. 


The  SVD  of  R  can  be  written  as 


R  =  UVUh  , 


where  U  is  a  unitary  matrix  and  V  is  a  diagonal  matrix.  Columns  of  U  are  eigenvectors  of  R ,  and  the 
elements  in  the  diagonal  of  V  are  the  corresponding  eigenvalues.  As  the  clutter  is  the  dominant 
component  in  the  received  signal,  the  eigenvectors  ul,u2,.-.,uM  corresponding  to  the  M  largest 
eigenvalues  can  be  reasonably  associated  to  the  clutter.  The  projection  of  the  received  signal  to  the 
orthogonal  subspace  of  the  clutter. 
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(22) 


M 

Sproj  -TUn,“n,  )Sc’ 

(=1 

results  in  clutter-suppressed  signal. 

In  (17),  because  the  number  of  neighboring  range  cells  is  usually  smaller  than  the  dimension  of  the 
variance  matrix  R  to  be  estimated,  it  is  rank  deficient.  By  considering  the  existence  of  thermal  noise,  the 
full-rank  covariance  matrix  Rt  of  clutter  and  noise  can  be  estimated  as 

tf,=/?  +  <72/,  (23) 

where  cr2  is  the  noise  variance  which  can  be  roughly  estimated,  and  /  is  the  identity  matrix.  Performing 
singular  value  decomposition  (SVD)  of  R]  yields 

/?,  =  t/.V.t/"  =  £  Ifruf1  =  £  X,P,  (24) 

1=1  1=1 

where  is  the  i-th  largest  eigenvalue  of  Rx ,  ui  is  its  eigen  vector,  and  P-  =  ufif  is  the  projection 
operator  to  the  subspace  generated  by  u( .  From  (21),  we  know  that  more  clutter  energy  distributed  in  the 
subspace  vector  ui  corresponding  to  a  larger  eigenvalue  Ai .  The  following  algorithm  provides  an 
efficient  way  to  remove  strong  clutter  without  any  knowledge  of  signal, 

SproJ-tfiWo  (25) 

1=1 

where  /(^)  is  a  weighting  function  that  takes  a  smaller  value  for  a  larger  value  of  Ai .  In  this  paper, 
f(%i)  is  chosen  as 

f(Ai)  =  yij.  (26) 

Therefore,  the  signal  vector,  after  the  adaptive  clutter  rejection,  becomes 


s„,i  =*fc  =(R+^/)-'se. 

*=l  i= 1  Ai 


The  noise  variance  estimate  <r2  controls  the  rejection  level  against  the  clutter  components. 
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5.  Simulation  Results 


In  this  section,  the  performance  of  the  proposed  algorithms  for  maneuvering  target  detection  is 
shown  by  some  simulation  results.  The  signal  data  coming  from  maneuvering  targets  is  generated  based 
on  the  signal  model  (1)  and  then  added  to  the  raw  OTHR  clutter  data.  The  radar  working  frequency  is 
20MHz.  There  are  54  range  cells  in  the  data.  The  coherent  integration  time  (CIT)  is  Tc  =  12.3  seconds. 

The  velocity  and  acceleration  of  targets  in  the  range  direction  are  from  40  m/s  and  3  m/s 2  respectively. 
The  signal  to  clutter  ratio  is  about  -53.5  dB. 

In  our  simulations,  the  following  steps  are  implemented.  For  the  received  signal,  matched  filtering 
and  range  compression  are  first  implemented  in  the  range  direction.  Then,  the  signal  subspace  clutter 
rejection  algorithm  is  applied  to  remove  the  clutter  where  y=l  is  used.  At  last  ACT  is  used  to  the  clutter- 
rejected  signal  for  target  detection. 

The  signal  waveforms  before  and  after  adaptive  signal  subspace  clutter  rejection  to  the  range  cell 
that  contains  target  are  shown  in  Figl.(a).  We  can  see  that  the  clutter  energy  is  removed  about  15  dB  by 
using  the  signal  subspace  algorithm.  The  results  can  also  be  verified  by  the  results  of  Fig.  1(b)  and 
Fig.  1(c),  which  are  the  mesh  of  range-Doppler  results  to  the  data  before  and  after  clutter  rejection  by 
adaptive  subspace  clutter  rejection  algorithm.  It  is  noted  that  the  clutter  suppression  at  edge  range  cells  is 
not  as  effective  as  the  other  cells  because  less  neighboring  range  cells  are  available  for  clutter  subspace 
estimation. 

The  processing  results  with  different  methods  are  shown  in  Fig.2.  The  target  can  not  be  detected 
from  the  range-Doppler  results  obtained  by  using  the  Fourier  transform  to  the  pre-clutter  rejection  data, 
which  is  shown  in  Fig2.(a).  Because  the  Fourier  transform  spreads  the  target  energy  of  the  maneuvering 
targets,  as  shown  in  Fig.  2(b),  the  target  is  still  undetectable  even  after  the  clutter  rejection  and  SCR 
enhancement.  Instead  of  the  Fourier  transform  in  Fig.2(b),  ACT  is  used  in  Fig.2(c).  The  target,  however, 
can  be  easily  detected  now  in  Fig.2(c)  at  range  2250  km  with  Doppler  about  -4  Hz.  The  amplitudes  of  the 
signal  in  the  range  cell  containing  the  target  are  shown  in  Figs.3(a)-(c)  for  the  Fourier  transform  without 
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clutter  rejection,  the  Fourier  transform  with  clutter  rejection,  and  the  ACT  with  clutter  rejection, 
respectively.  From  Fig.3(a),  we  can  see  that  the  clutter  amplitude  of  the  main  lobe  around  0  Hz  is  about 
20  dB  higher  than  that  of  the  side  lobes  including  the  region  around  -5  Hz  where  the  target  is  located.  In 
Fig.3(b),  although  the  clutter  amplitude  of  the  main  lobe  is  reduced  by  about  15  dB  and  the  side  lode  is 
reduced  about  5  to  10  dB,  the  target  still  can  not  be  detected.  But  in  Fig.3(c),  the  target  energy  is  focused. 
The  amplitude  of  the  target  signal  is  about  4  dB  higher  than  that  of  the  clutter  in  target’s  neighboring 
frequency  bands. 


6.  Conclusion 

In  this  paper,  an  adaptive  clutter  rejection  algorithm  was  proposed  to  maneuvering  target  detection 
in  OTHR  systems.  This  algorithm  can  reduce  clutter  energy  by  about  15  to  20  dB  with  negligible 
distortion  to  the  waveform  of  the  signal  returned  from  maneuvering  targets.  An  adaptive  chirplet 
transform  algorithm  was  applied  to  the  clutter-mitigated  signal  for  improved  Doppler  processing. 
Simulation  results  showed  that  the  proposed  method  substantially  enhances  the  target  detection  ability. 
Particularly,  several  simulation  examples  showed  that  the  proposed  method  can  successfully  detect  weak 
target  signals  where  other  methods  can  not  be  used  directly  with  adaptive  chirplet  transform  for 
maneuvering  target  detection. 
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Abstract 

Bilinear  synthesis  of  nonstationary  signals  impinging  on  a  multi-antenna  receiver  has  been  recently 
introduced.  The  distinction  in  the  spatial  signatures  of  the  sources  provides  a  vehicle  to  reduce  noise 
and  source  signal  interactions  in  the  time-frequency  domain,  and  hence  improves  signal  synthesis.  In 
this  letter,  we  utilize  another  form  of  diversity  for  enhanced  source  time-frequency  signal  representations. 
It  is  shown  that  cross-polarization  antennas  can  be  used  to  mitigate  crossterms  via  simple  polarization 
averaging. 
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I.  Introduction 


Time-frequency  distributions  (TFDs)  have  been  found  useful  in  the  analysis  and  classi¬ 
fication  of  nonstationary  signals  [1],  [2].  In  [3],  it  is  shown  that  the  array  manifold  can  be 
used  to  improve  syntheses  of  signals  with  rapid  time- varying  frequency  characteristics.  In 
essence,  averaging  TFDs  across  different  array  sensors  trades  off  the  spatial  dimension  for 
enhanced  auto-source  TFDs.  Spatial  averaging  mitigates  the  cross-source  time-frequency 
(t-f)  terms  as  well  as  reduces  the  noise  contribution. 

When  the  receiver  is  not  equipped  with  an  antenna  array,  or  the  array  is  of  small  aper¬ 
ture,  the  spatial  averaging  of  TFDs  proposed  in  [3]  will  no  longer  be  effective  or  applicable. 
A  possible  alternative  is  to  use  cross-polarization  antennas  where  the  polarization  dimen¬ 
sion  can  be  utilized  to  enhance  t-f  signature  estimation  and  subsequently  leads  to  improved 
signal  synthesis  performance.  The  polarization-based  t-f  signal  synthesis  can  be  used  for  a 
single  as  well  as  multiple  antennas.  In  this  letter,  we  restrict  our  discussion  to  the  simple 
case  of  a  single  pair  of  cross-polarization  antennas.  The  generalization  to  applications  of 
multi-sensor  receivers  is  straightforward  and  is  addressed  in  [4]. 

Signal  polarization  properties  have  been  commonly  utilized  in  wireless  communications 
and  synthetic  aperture  radars  (SARs)  [5],  [6].  Distinct  polarization  signatures  of  differ¬ 
ent  sources  can  be  observed  when  the  sources  have  different  transmitter  polarizations  or 
distinct  channel  characteristics. 

This  letter  is  organized  as  follows.  The  signal  model  is  presented  in  Section  II.  Section 
III  proposes  the  polarization  averaging  for  t-f  signal  synthesis.  The  analogy  between  the 
proposed  method  and  the  array  averaging  is  also  considered.  Numerical  simulations  are 
given  in  Section  IV. 

II.  Signal  Model 

The  discrete-time  data  received  at  a  cross-polarization  antenna,  which  receives  two 
orthogonal  polarizations  (e.g.,  vertical  and  horizontal  polarizations),  is  expressed  in  the 
following  vector  format 

x(f)  =  [zw(t)  a:f,?](t)]T,  (1) 

where  M  and  M  represent  the  two  orthogonal  polarizations,  and  T  denotes  transpose. 
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The  following  expression 

OO  OO 

Dx[i)x[k](t,  f)  =  j  J  <j)(t-  u,  T)x[l](t  +  ^)(x[k](t  -  ^))*e~32vftdudr  (2) 

—  OO  — OO 

defines  the  auto-  ( i  =  k )  and  cross-polarized  (i  ^  k)  TFDs  of  the  two  polarizations, 
where  t  and  /  are  the  time  and  the  frequency  indeces,  respectively,  and  r)  is  the 
time-frequency  kernel  [7].  Each  of  i  and  k  takes  either  value  of  the  polarization  index  p 
or  q.  The  auto-  and  cross-polarized  TFDs  can  be  combined  to  form  the  following  2x2 
polarization  TFD  matrix, 

OO  OO 

Dxx(f,  f)=  j  f  4>(t-  u,  r)x(t  +  ^)xH(Z  -  ^)e~j2nftdudT.  (3) 

—  OO  — OO 

where  superscript  H  denotes  transpose  conjugation.  The  diagonal  entries  of  Dxx(t,  /)  are 
the  auto-polarized  TFDs,  whereas  the  off-diagonal  elements  are  cross-polarized  TFDs. 

Assume  L  source  signals  =  1 are  incident  on  the  antenna.  The  received 

data  for  each  polarization  is  the  linear  combination  of  the  same  polarization  components 
of  the  source  signals  and  noise.  That  is, 

xW(t)  =  ofW*)  +  *W(*)»  i  =  P,Q,  (4) 

/=1 

where  aft  represents  the  mixing  coefficient  of  the  Zth  source  along  the  ith  polarization, 
and  n^(t)  is  the  noise  component  at  the  same  polarization.  In  the  vector  form,  x(t)  can 
be  decomposed  into  the  following  terms 

L 

x(i)  =  y  (t)  +  n(t)  =  As  (t)  +  n(t)  =  ^  a ist(t)  +  n  (t),  (5) 

i=i 

where sq  =  [o^1  a[9l]T,  A  =  [ai, . . . ,  aL],  s(f)  =  [si(t), . . .  ,Si(t)]r,  and  n(t)  =  [nw(t)  nM(t)]r. 
Because  of  the  ambiguity  with  respect  to  the  signal  strength  and  the  propagation  atten¬ 
uation,  it  is  convenient  to  assume  that  Ha/JH  =  2,  l  =  1, •••,£,  and  the  propagation 
attenuation  scalar  is  absorbed  in  si(t).  The  noise  elements  are  modelled  as  stationary  and 
white  complex  Gaussian  processes  with  zero  mean  and  variance  <r2  in  each  polarzation, 
i.e., 

E  [n(t  +  r)nH(l)]  =  a2^(r)I2,  (6) 

where  S(r)  is  the  Kronecker  delta  and  IN  denotes  the  N  x  N  identity  matrix. 
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III.  Polarization  Averaging 

It  is  clear  from  Section  II  that  the  signal  model  of  the  cross-polarization  antenna  case  is 
similar  to  that  of  a  two-antenna  array.  The  only  difference  is  that  the  source  polarization 
vector  a*  is  used  in  place  of  the  source  spatial  signature,  or  steering  vector.  Accordingly, 
polarization  averaging  can  be  equally  effective  as  spatial  averaging  in  mitigating  the  TFD 
crossterms. 

From  eqns.  (2)  and  (4),  the  auto-polarization  TFD  of  (t)  is  given  by 

L  L 

Dx[i]x[i](t,  /)  =  EE  4](alS)*DSlSm(t,  f )  +  Dn[i]nli](t,  /),  i  =  p,q,  (7) 

1= 1  m—  1 

where  DSlSjn(t,f )  represents  the  auto-source  TFD  (if  l  =  m)  or  the  cross-source  TFD  (if 
l  y £  m).  The  presence  of  cross-source  terms  often  obscures  the  true  power  localization  over 
time  and  frequency. 

Averaging  the  auto-polarization  TFDs  over  the  two  polarization  branches  yields 

w  (t,f)  = 

^  i=p 

=  EE  (IT,  af  (am)*)  DSlSm  (t,  /)  +  i  E  PnMnM  (*,  /)] 

i=lj=l\*i=p  /  Z  i=P 

L  L  / 1  \  i  Q 


E E  (h&«)  /)  +  i E /)]  • 

1=1  7=1  L  i=p 


In  eqn.  (8),  a^a*  is  the  inner  product  of  the  polarization  signatures  a m  and  a i.  Define  the 
polarization  correlation  coefficient 

Am  =  2&maO  (9) 

Accordingly,  eqn.  (8)  can  be  expressed  as 

W(i.  /)  =  E  E  /)  +  \  E  /)] .  (10) 

i=l  j= 1  Z  i=p 

The  above  equation  shows  that  W (t,  /)  is  a  linear  combination  of  the  auto-  and  cross¬ 
polarization  TFDs  of  all  signal  arrivals.  It  is  straightforward  to  show  that  for  the  Ith  and 
the  mth  sources, 


I  Am  |  <  1,  if  l^m  and  Am  =  1,  if  l  =  m, 
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indicating  that  the  constant  coefficients  in  (10)  for  the  auto-polarization  TFDs  are  always 
greater  than,  or  at  least  equal  to,  those  for  the  cross-polarization  TFDs.  For  sources  with 
distinct  polarizations,  \Pim\  <C  1,  leading  to  significant  suppression  of  crossterms,  and 
thereby  enhancing  the  signal  signature  estimation. 

An  interesting  case  arises  when  two  signals  have  orthogonal  polarization  signatures,  i.e., 
(3im  =  0  for  l  ±  m.  In  this  case,  the  crossterms  between  these  two  source  signals  will  be 
entirely  eliminated  and  only  the  autoterms  will  be  maintained. 

The  t-f  kernel  in  eqns.  (2)  and  (3),  which  introduces  temporal  averaging  of  the  local 
autocorrelation  functions  at  consecutive  time  samples,  can  be  selected  to  reduce  the  TFD 
noise  effect  for  the  single  antenna  case,  as  discussed  in  [8],  [9].  However,  even  without 
kernel  smoothing,  the  polarization  averaging  in  eqn.  (10),  similar  to  spatial  averaging 
[10],  decreases  the  noise  variance  and  its  interaction  with  the  signal  components  beyond 
that  achieved  in  a  single  antenna  (polarization)  case.  Once  the  polarization  averaging 
is  performed  and  the  t-f  signature  is  identified,  we  can  then  proceed  with  the  bilinear 
syntheses  using  the  methods  described  in  [2]. 

It  is  noted  that,  although  the  model  used  allows  for  L  source  signals  to  be  present,  there 
are  only  two  dimensions  of  polarization  diversity  for  a  single  cross-polarization  antenna. 
Therefore,  when  L  >  2,  while  the  crossterms  between  different  source  signals  can  still  be 
substantially  mitigated,  it  becomes  impossible  to  completely  eliminate  all  the  crossterms 
unless  more  sensors  are  used. 


IV.  Simulation  Results 


In  this  section,  we  provide  computer  simulations  to  demonstrate  the  improvement  gained 
by  the  proposed  technique  in  the  reduction  or  elimination  of  cross  terms  and  signal  syn¬ 
thesis.  Two  high-order  frequency  modulated  signals  are  considered  on  a  dual-polarization 
dipole.  Their  polarizations  are  assumed  to  be  orthogonal,  with  the  following  mixing  ma- 


The  length  of  the  signal  sequence  is  set  to  N  —  256.  The  additive  noise  is  zero  mean, 
Gaussian  distributed,  and  white.  The  input  signal-to-noise  ratio  (SNR)  is  3  dB. 
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With  the  presence  of  high-level  noise  and  close  t-f  signatures,  it  is  very  difficult  to 
identify  these  t-f  signatures  when  only  a  single-polarization  sensor  is  used.  Fig.  1  shows 
the  extended  discrete-time  Wigner-Ville  distribution  (EDTWVD)  [11]  of  the  data  received 
at  the  vertically  polarized  antenna.  However,  as  evident  from  Fig.  2,  the  t-f  signatures 
of  the  two  signals  can  be  revealed  when  polarization  averaging  is  applied.  Fig.  2  shows 
that  the  crossterm  between  the  two  signals  are  completely  eliminated  and  the  variance  of 
noise  terms  is  reduced.  Masking  the  first  signal  and  applying  standard  signal  synthesis 
techniques  yield  a  high  quality  signal  recovery.  Fig.  3  shows  the  TFD  of  the  synthesized 
signal  waveform  of  the  first  signal. 


TFD  at  a  sensor 


time 


Fig.  1.  EDTWVD  computed  from  the  signal  received  at  the  vertical  polarization 

antenna. 


averaged  TFD 


time 


Fig.  2.  EDTWVD  averaged  over  two  polarizations. 
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TFD  of  synthesized  signal 


Fig.  3.  EDTWVD  of  the  synthesized  waveform  of  the  first  signal. 


V.  Conclusion 

Polarization  averaging  allows  effective  crossterm  reduction  and  autoterm  enhancement, 
aiding  source  time-frequency  signature  estimations  and  waveform  recovery.  Averaging 
TFDs  across  polarizations  can  be  performed  concurrently  with  TFD  averaging  across  the 
array,  thereby  utilizing  both  spatial  and  polarization  diversity  in  syntheses  of  nonstationary 
signals.  However,  polarization  averaging  can  be  applied  alone  if  the  difference  in  the  source 
spatial  signatures  is  insufficient  for  crossterm  reduction,  or  the  receiver  is  not  equipped 
with  antenna  arrays. 
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Abstract 

This  chapter  presents  a  comprehensive  treatment  of  the  hybrid  area  of  time-frequency  distribution 
(TFD)  and  array  signal  processing.  The  application  of  quadratic  time-frequency  distributions  to  sensor 
signal  processing  has  recently  become  of  interest,  and  it  was  necessitated  by  the  need  to  address  important 
problems  related  to  processing  nonstationary  signals  incident  on  multi-antenna  receivers.  Over  the  past 
few  years,  major  contributions  have  been  made  to  improve  direction  finding  and  blind  source  separation 
using  time-frequency  signatures.  This  improvement  has  cast  quadratic  time-frequency  distributions  as 
a  key  tool  for  source  localization  and  signal  recovery,  and  put  bilinear  transforms  at  equal  footing  with 
second-order  and  higher-order  statistics  as  bases  for  effective  spatial-temporal  signal  processing.  This 
chapter  discusses  the  advances  made  through  time-frequency  analysis  in  direction-of-arrival  estimation, 
signal  synthesis,  and  near-field  source  characterization. 
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I.  Introduction 


Time-frequency  distributions  (TFDs)  are  used  in  various  applications,  including  speech, 
biomedicine,  sonar,  radar,  GPS,  and  geophysics.  Over  the  past  two  decades,  most  of  the 
work  on  quadratic  TFDs  has  focused  on  the  mono-component  and  multi-component  tem¬ 
poral  signal  structures  and  the  corresponding  time-frequency  (t-f)  signatures.  This  work 
has  led  to  major  advances  in  nonstationary  signal  analysis  and  processing.  Information 
on  the  signal  instantaneous  frequency  and  instantaneous  bandwidth  obtained  from  the  t-f 
domain  has  allowed  improved  separation,  suppression,  classification,  identification,  and 
detection  of  signals  with  time  varying  spectra  [1],  [2],  [3],  [4],  [5],  [6]. 

Applications  of  the  quadratic  distributions  to  array  signal  processing  began  to  flourish 
in  the  mid-nineties.  The  main  objective  was  to  enhance  direction-finding  and  blind  source 
separation  of  nonstationary  signals  using  their  t-f  signatures.  Another  important  but 
different  objective  was  to  characterize  near-field  and  far-field  emitters  based  on  their  spatial 
signatures.  In  order  to  achieve  both  objectives,  new  definitions  and  generalization  of 
quadratic  distributions  were  in  order. 

The  spatial  time-frequency  distribution  (STFD)  has  been  introduced  to  describe  the 
mixing  of  nonstatioanry  signals  at  the  different  sensors  of  the  array.  The  relationship 
between  the  TFDs  of  the  sensors  to  the  TFDs  of  the  individual  source  waveforms  is  defined 
by  the  steering,  or  the  array,  matrix,  and  was  found  to  be  similar  to  that  encountered  in 
the  traditional  data  covariance  matrix  approach  to  array  processing.  This  similarity  has 
allowed  a  rapid  progress  in  nonstationary  array  processing  from  the  TFD  perspective  [7]. 

This  chapter  discusses  two  fundamental  formulations  to  incorporate  the  spatial  informa¬ 
tion  into  quadratic  distributions.  One  formulation  is  based  on  STFDs  and  the  localization 
of  the  signal  arrivals  in  the  t-f  domain.  The  corresponding  analysis,  theory,  and  applica¬ 
tions  are  covered  in  the  first  six  sections  of  this  Chapter.  Section  VII  deals  with  another 
formulation,  in  which  the  quadratic  distribution  of  the  spatial  signal  across  the  array  is 
computed.  This  sensor-angle  distribution  (SAD)  localizes  the  source  angle  at  each  sensor, 
and  is  the  dual  in  sensor  number  and  angle  to  Cohen’s  class  of  time-frequency  distributions. 
The  SAD  is  particularly  appropriate  for  characterizing  sources  and  scatter  in  the  near-field 
of  an  array.  Sources  arriving  from  the  far-field  have  the  same  angle  at  each  sensor.  In  con- 
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trast,  sources  in  the  near-field  have  differing  angle  at  each  sensor,  and  a  full  sensor-angle 
distribution  provides  a  complete  characterization  of  the  near-field  source  and  scatter  en¬ 
vironment.  Knowledge  of  this  characterization  can  be  important  when  calibrating  arrays 
of  sensors  placed  in  a  non-homogeneous  environment,  such  as  a  radar  or  communications 
array  deployed  in  a  built-up  environment  and  surrounded  by  other  metallic  structures  that 
are  not  part  of  the  array. 

II.  Spatial  Time-Frequency  Distributions 

A.  Signal  Model 

In  narrowband  array  processing,  when  n  signals  arrive  at  an  m-element  (sensor)  array, 
the  linear  data  model 

x(£)  =  y(t)  +  n(t)  =  Ad(£)  +  n(t)  (1) 

is  commonly  assumed,  where  the  m  x  n  spatial  matrix  A  =  [ai,  •  •  • ,  a„]  represents  the 
mixing  matrix  or  the  steering  matrix.  In  direction  finding  problems,  we  require  A  to  have 
a  known  structure,  and  each  column  of  A  corresponds  to  a  single  arrival  and  carries  a 
clear  bearing.  For  blind  source  separation  problems,  A  is  a  mixture  of  several  steering 
vectors,  due  to  multipaths,  and  its  columns  may  assume  any  structure. 

The  mixture  of  the  signals  at  each  sensor  renders  the  elements  of  the  mxl  data  vector 
x(t)  to  be  multicomponent  signals,  whereas  each  source  signal  cf*(£)  of  the  n  x  1  signal 
vector  d(£)  is  typically  a  monocomponent  signal.  n(t)  is  an  additive  noise  vector  whose 
elements  are  modeled  as  stationary,  spatially  and  temporally  white,  zero-mean  complex 
Gaussian  random  processes,  independent  of  the  source  signals.  That  is, 

E[n(t  +  r)nH(£)]  =  aS(r) I  and  E[n(t  +  r)nT(f)]  =  0  for  any  r  (2) 

where  S(t)  is  the  delta  function,  I  denotes  the  identity  matrix,  a  is  the  noise  power  at 
each  sensor,  superscript  H  and  T,  respectively,  denote  conjugate  transpose  and  transpose, 
and  E(-)  is  the  statistical  expectation  operator. 

B.  Spatial  Time- Frequency  Distributions 

We  first  review  the  definition  and  basic  properties  of  the  spatial  time-frequency  distribu¬ 
tions  (STFDs).  STFDs  based  on  Cohen’s  class  of  time-frequency  distribution  (TFD)  were 
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introduced  in  [8]  and  their  applications  to  direction  finding  and  blind  source  separation 
have  been  discussed  in  [9],  [10]  and  [8],  [11],  respectively. 

The  discrete  form  of  TFD  of  a  signal  x(t)  is  given  by 

OO  00 

Dxx(t,  f)=  Yj  E  l)x(f  +  v  +  Ox*(*  +  v  -  (3) 

v=-oo  l=~ OO 

where  (f>(v,  l )  is  a  kernel  function  and  *  denotes  complex  conjugate.  The  STFD  matrix  is 
obtained  by  replacing  x(t)  by  the  data  snapshot  vector  x(t), 

OO  OO 

Dxx(f,  /)  =  Yj  E  Ox(*  +  v  +  +  v  -  /)e_j47r/f,  (4) 

v=~oo/=— OO 

Substitute  (1)  into  (4),  we  obtain 

Dxx(t,  /)  =  Dyy(t,  /)  +  Dyn(i,  /)  +  Dny(f,  /)  +  Dnn(t,  /).  (5) 

We  note  that  Dxx(t, /),  Dyy (<,/),  Dyn(t, /),  Dny(i,/),  and  Dnn(t,/)  are  matrices  of 
dimension  mxm.  Under  the  uncorrelated  signal  and  noise  assumption  and  the  zero-mean 
noise  property,  the  expectation  of  the  crossterm  STFD  matrices  between  the  signal  and 
noise  vectors  is  zero,  i.e.,  E  [Dyn(t,  /)]  =  E  [Dny(i,  /)]  =  0.  Accordingly, 

E  [Dxx(t,  /)]  =  Dyy(t,  f)  +  E  [Dnn(t,  /)] 

=  ADdd(t,  f)AH  +  E  [Dnn(t,  /)] ,  (6) 

where  the  source  TFD  matrix 

OO  OO 

Ddd (*,  /)  =  E  E  Hv,  l)d(t  +  v  +  l)dH(t  +  v-  l)e-^fl  (7) 

V=-00  l=- OO 

is  of  dimension  n  x  n.  For  narrowband  array  signal  processing  applications,  the  mixing 
matrix  A  holds  the  spatial  information  and  maps  the  auto-  and  cross-TFDs  of  the  source 
signals  into  auto-  and  cross-TFDs  of  the  data. 

Equation  (6)  is  similar  to  the  formula  that  has  been  commonly  used  in  direction  finding 
and  blind  source  separation  problems,  relating  the  signal  correlation  matrix  to  the  data 
spatial  correlation  matrix.  In  the  above  formulation,  however,  the  correlation  matrices 
are  replaced  by  the  STFD  matrices.  The  well  established  results  in  conventional  array 
signal  processing  could,  therefore,  be  utilized  and  key  problems  in  various  applications  of 
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array  processing,  specifically  those  dealing  with  nonstationary  signal  environments,  can 
be  approached  using  bilinear  transformations. 

Initially,  only  the  time-frequency  (t-f)  points  in  the  autoterm  regions  of  TFD  are  con¬ 
sidered  for  STFD  matrix  construction.  The  autoterm  region  refers  to  the  t-f  points  along 
the  true  instantaneous  frequency  (IF)  of  each  signal.  The  crossterms,  which  can  intrude 
on  the  autoterms  through  the  power  in  their  mainlobes  or/and  sidelobes,  were  avoided. 
This  intrusion  depends  on  the  signal  temporal  structures  and  the  window  size.  Recently, 
the  crossterms  have  also  been  utilized  and  integrated  into  STFDs.  The  effect  of  crossterms 
on  direction  finding  and  blind  source  separation  will  be  discussed  in  IV-C  and  V-C,  re¬ 
spectively.  In  the  other  parts  of  this  chapter,  it  is  assumed  that  the  t-f  points  reside  in  an 
autoterm  region,  which  has  negligible  crossterm  effect. 

C.  Joint- Diagonalization  and  Time- Frequency  Averaging 

In  the  rest  of  this  chapter,  we  will  address  the  application  of  STFDs  to  direction  finding 
and  blind  source  separation.  These  applications  are  based  on  the  eigendecomposition  of 
the  STFD  matrix.  In  direction  finding,  the  source  TFD  matrix  must  be  full  rank  (Section 
V),  whereas,  to  perform  blind  source  separation,  the  source  TFD  matrix  must  be  diagonal 
(Section  IV).  For  either  case,  the  STFD  matrix  of  the  data  vector  should  be  full  column 
rank. 

It  is  noted  that  the  relationship  (6)  holds  true  for  every  (t,  f)  point.  In  order  to  en¬ 
sure  the  full  column  rank  property  of  the  STFD  matrix  as  well  as  to  reduce  the  effect  of 
noise,  we  consider  multiple  t-f  points,  instead  of  a  single  one.  This  allows  more  informa¬ 
tion  of  the  source  signal  t-f  signatures  to  be  included  into  their  respective  eigenstructure 
formulation,  and  as  such  enhances  direction  finding  and  source  separation  performance. 
Joint-diagonalization  and  t-f  averaging  are  the  two  main  approaches  that  have  been  used 
for  this  purpose  [8],  [9],  [12]. 
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C.l  Joint  Diagonalization 

The  JD  can  be  explained  by  first  noting  that  the  problem  of  the  diagonalization  of  a 
single  n  x  n  normal  matrix  M  is  equivalent  to  the  minimization  of  the  criterion  [13] 

C'(M,V)  =  -^|vfMvi|2  (8) 

i 

over  the  set  of  unitary  matrices  V  =  [vi,---,v„].  Hence,  the  JD  of  a  set  {Mfc|A;  = 
1,  •  •  • ,  K}  of  K  arbitrary  n  x  n  matrices  is  defined  as  the  minimization  of  the  following 
JD  criterion: 

C(V)  =  -£C(M,V)  =  -££|v'fMivi|2  (9) 

k  k  i 

under  the  same  unitary  constraint.  It  is  important  to  note  that  the  above  definition  of 
JD  does  not  require  the  matrix  set  under  consideration  to  be  exactly  and  simultaneously 
diagonalized  by  a  single  unitary  matrix.  This  is  because  we  do  not  require  the  off-diagonal 
elements  of  all  the  matrices  to  be  cancelled  by  a  unitary  transform;  a  joint  diagonalizer 
is  simply  a  minimizer  of  the  criterion.  If  the  matrices  in  M  are  not  exactly  joint  diag- 
onalizable,  the  criterion  cannot  be  zeroed,  and  the  matrices  can  only  be  approximately 
joint  diagonalized.  Hence,  an  (approximate)  joint  diagonalizer  defines  a  kind  of  average 
eigenstructure. 

C.2  Joint  Block-Diagonalization 

For  direction  finding  methods  such  as  t-f  MUSIC,  the  source  TFD  matrix  should  not  be 
singular  but  not  necessarily  diagonal.  In  this  case,  the  joint  block-diagonalization  ( JBD) 
is  used  to  incorporate  multiple  t-f  points  rather  than  JD  [9].  The  JBD  is  achieved  by  the 
maximization  under  unitary  transform  of  the  following  criterion 

C(U)  =  ^£|ufMuJ2  (10) 

k  i,l 

over  the  set  of  unitary  matrices  U  =  [ui,  •  •  • ,  un]. 

C.3  Time-Frequency  Averaging 

Time-frequency  averaging  is  a  linear  operation  that  adds  the  STFDs  over  a  t-f  region 
where,  typically,  the  desired  signal  is  highly  localized  and  the  crossterms  are  negligible. 
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The  averaged  STFD  is  defined  as 


0  =  7  E  D „(*,/).  (11) 

where  Q  is  the  t-f  region  of  interest  and  A  is  a  normalization  constant  which,  for  example, 
can  be  chosen  as  the  total  number  of  (t,  f)  points  in  the  region  £1  The  eigendecomposition 
of  D  is  addressed  in  III-C. 


III.  Properties  of  STFDs 

To  understand  the  properties  of  STFDs,  we  consider  the  case  of  frequency  modulated 
(FM)  signals  and  the  simplest  form  of  TFD,  namely,  the  pseudo  Wigner-Ville  distribution 
(PWVD)  [14].  The  consideration  of  FM  signals  is  motivated  by  the  fact  that  these  signals 
are  uniquely  characterized  by  their  IFs,  and  therefore,  they  have  clear  t-f  signatures  that 
can  be  utilized  by  the  STFD  approach.  Also,  FM  signals  have  constant  amplitudes. 

The  FM  signals  can  be  modeled  as 

d(t)  =  [di(t), ...,  dn(t)]T  =  [£>ie^l(t), ...,  Dnejxpn{t)}T ,  (1) 

where  Di  and  ipi(t)  are  the  fixed  amplitude  and  time- varying  phase  of  ith  source  signal. 
For  each  sampling  time  t ,  dt(t)  has  an  instantaneous  frequency  /*(f)  =  dipi(t) /(2ndt). 

The  discrete  form  of  PWVD  of  a  signal  x(t),  using  a  rectangular  window  of  odd  length 
L,  is  a  special  case  of  (3)  and  is  given  by 

(£-l)/2 

Dxx(t,  f)=  x(t  +  t )x* (t  -  T)e~jA*fT.  (2) 

r=-(L-l)/2 

Similarly,  the  spatial  pseudo  Wigner-Ville  distribution  (SPWVD)  matrix  is  obtained  by 
replacing  x(t)  by  the  data  snapshot  vector  x(£), 

(L-l)/2 

D xx(t,f)=  ^2  x(£  +  r)xH(t  -  r)e_i47r/T.  (3) 

r=-(L-l)/2 

A.  Subspace  Analysis  for  FM  Signals 

Analysis  of  the  eigendecomposition  of  the  STFD  matrix  is  closely  related  to  the  anal¬ 
ysis  of  subspace  decomposition  of  the  covariance  matrix  [15].  Before  elaborating  on  this 
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relationship,  we  present  the  case  of  FM  signals  using  the  conventional  covariance  matrix 
approach. 

In  equation  (1),  it  is  assumed  that  the  number  of  sensors  is  greater  than  the  number  of 
sources,  i.e.,  m  >  n.  Further,  matrix  A  is  full  column  rank.  We  further  assume  that  the 
correlation  matrix 

Rxx  =  £[x(f)xff(t)]  (4) 

is  nonsingular,  and  the  observation  period  consists  of  N  snapshots  with  N  >  m.  Under 
the  above  assumptions,  the  correlation  matrix  is  given  by 

Rxx  =  £[x(i)x*(i)]  =  ARddAH  +  ffl,  (5) 

where  Rdd  =  E[d(t)dH(t)}  is  the  source  correlation  matrix. 

Let  Ax  >  A2  >  •  •  •  >  An  >  A„+1  =  A„+2  =  •  •  •  =  Am  =  a  denote  the  eigenvalues  of  Rxx. 
It  is  assumed  that  A*,  i  =  1,  •  •  • ,  n,  are  distinct.  The  unit-norm  eigenvectors  associated 
with  Ai,...,An  constitute  the  columns  of  matrix  S  =  [si,...,s„],  and  those  corresponding 
to  An+i, ...,  Am  make  up  matrix  G  =  [gi, ...,  gm_„j.  Since  the  columns  of  A  and  S  span  the 
same  subspace,  then  AHG  =  0. 

In  practice,  Rxx  is  unknown,  and  therefore  should  be  estimated  from  the  available  data 
samples  (snapshots)  x(i),  i  =  1,2, ...,  N.  The  estimated  correlation  matrix  is  given  by 

=  -^r]Tx(i)xH(i).  (6) 

ly  i=l 

Let  {§i, sn,  gi, ...,  grn— n}  denote  the  unit-norm  eigenvectors  of  Rxx,  arranged  in  the 
descending  order  of  the  associated  eigenvalues,  and  let  S  and  G  denote  the  matrices 
defined  by  the  set  of  vectors  {§,}  and  {gi},  respectively.  The  statistical  properties  of 
the  eigenvectors  of  the  sample  covariance  matrix  Rxx  for  signals  modeled  as  independent 
processes  with  additive  white  noise  are  given  in  [15]. 

We  assume  that  the  transmitted  signals  propagate  in  a  stationary  environment  and 
are  mutually  uncorrelated  over  the  observation  period  1  <  t  <  N.  Subsequently,  the 
corresponding  covariance  matrices  are  time-independent.  Under  these  assumptions, 

=  0  for*#J,  i,i  =  l,-,n.  (7) 

iV  k=l 
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In  this  case,  the  signal  correlation  matrix  is 

1  T 

Rdd  =  lim  -  ^2  d(t)dH(t)  =  diag  [Df,  i  =  1,2, ...,  n]  (8) 

where  diag[ ■]  is  the  diagonal  matrix  formed  with  the  elements  of  its  vector  valued  argu¬ 
ments.  From  the  above  assumptions,  we  have  the  following  Lemma. 

Lemma  1  [14]:  For  uncorrelated  FM  signals  with  additive  white  Gaussian  noise, 
a)  The  estimation  errors  (s2  —  Sj)  are  asymptotically  (for  large  N )  jointly  Gaussian 
distributed  with  zero  means  and  covariance  matrices  given  by 


E  [(Si  -  s i){Sj  -  s,-)*] 


_  n  \  ,  \  —  m—n  \ 

O  -f  Afc  V  jj  (  Ai  H  x 

N  h  (A*  -  A,)2  SfcSfc  ti  (<r  -  Aj)2gfcgfc  ^ 


E  [(Sf  -  Si)(8j  -  Sj)T ] 


J>1  <T  (Aj  +  A  j  (7 )  rjp 


N  (A,  -  \i) 


fVF(l  -  iy). 


where 


d«  =  {'’  '  J' 

l  0,  j. 

b)  The  orthogonal  projections  of  {g,}  onto  the  column  space  of  S  are  asymptotically 
(for  large  N)  jointly  Gaussian  distributed  with  zero  means  and  covariance  matrices  given 

by 


E  [(ss*gi)  (SS*g $ 


hi  a 


H  f  dcf  *  TTf 


N  Lti  (*  -  A*)2 


E  [(SSffgj)  (SSffgi)T]  =  0  for  all  i,  j.  (12) 

Equations  (9)  and  (10)  hold  strong  similarity  to  those  of  [15].  The  only  difference  is  that  the 
term  (AjAfc)  in  [15]  is  replaced  by  cr(A;  +  A*  —  a)  in  (9)  and  (10),  due  to  the  uncorrelation 
property  (7).  Equations  (11)  and  (12)  are  identical  to  those  derived  in  reference  [15]. 

B.  SNR  Enhancement 

The  TFD  maps  one-dimensional  signals  in  the  time  domain  into  two-dimensional  signals 
in  the  t-f  domain.  The  TFD  property  of  concentrating  the  input  signal  around  its  IF 
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while  spreading  the  noise  over  the  entire  t-f  domain  increases  the  effective  SNR  and  proves 
valuable  in  the  underlying  problem. 

The  ith  diagonal  element  of  PWVD  matrix  Ddd  (t,  /)  is  given  by 

(L-l)/2 

Ddidi(tJ)=  E  DzertMt+T)-Mt-T)}-j4«fT  (13) 

t=-(L- l)/2 

Assume  that  the  third-order  derivative  of  the  phase  is  negligible  over  the  window  length 
L,  then  along  the  true  t-f  points  of  the  zth  signal,  fi(t)  =  dipi(t)/(2ndt),  and  ^(t  +  r)  — 
ipi (t  —  r)  —  47t fi{t)r  =  0.  Accordingly,  for  ( L  —  l)/2  <  t  <  N  —  ( L  —  l)/2, 

(L- 1)/2 

Ddidi(t,  =  E  Di  =  LDl  (14) 

r=-(L- 1)/2 

Similarly,  the  noise  SPWVD  matrix  Dnn(t,  /)  is 

(i-i)/2 

Dnn (*,/)=  E  n(t  +  r)nH (t  -  r)e~j4vfT .  (15) 

t=-(L- l)/2 

Under  the  spatially  and  temporally  white  assumptions,  the  statistical  expectation  of 
Dnn(t,  /)  is  given  by 

(L-l)/2 

£(D„(t./)]=  2  £[n(<  +  r)nH(J-T)]e--’4^=rt  (16) 

r=-(i-l)/2 

Therefore,  when  we  select  the  t-f  points  along  the  t-f  signature  or  the  IF  of  the  ith  FM 
signal,  the  SNR  in  the  model  (6)  becomes  LDf/a,  which  has  an  improved  factor  L  over 
the  one  associated  with  model  (5).  The  IF  of  the  FM  signals  can  be  estimated  from  the 
employed  TFD,  or  using  any  appropriate  IF  estimator.  It  is  noted,  however,  that  the 
STFD  equation  (6)  provides  a  natural  platform  for  the  direct  incorporation  of  any  a  priori 
information  or  estimates  of  the  IF  into  direction-of-arrival  (DOA)  estimation. 

The  PWVD  of  each  FM  source  has  a  constant  value  over  the  observation  period,  pro¬ 


viding  that  we  leave  out  the  rising  and  falling  power  distributions  at  both  ends  of  the 
data  record.  For  convenience  of  analysis,  we  select  those  N'  —  N  —  L  -{■  1  t-f  points  of 
constant  distribution  value  for  each  source  signal.  In  the  case  where  the  STFD  matrices 
are  averaged  over  the  t-f  signatures  of  na  sources,  i.e.,  a  total  of  n0N'  t-f  points,  the  result 


is  given  by 


D  = 


1 

n0N' 


n0  Nr 

EE  Dxx(*t,/9,i(*i)), 


q=  1  i—1 


(17) 
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where  fq,i(U )  is  the  IF  of  the  qth  signal  at  the  ith  time  sample.  x(f)  is  an  instantaneous 
mixture  of  the  FM  signals  rfj(t),  i  =  1,  •  •  • ,  n,  hence  features  the  same  IFs.  The  expectation 
of  the  averaged  STFD  matrix  is 

1  n0  N' 

D  =  £[6]  =  -lf7EEBtDx,(<i,  /„,<((.))! 

I  n0  T 

=  -E[iflVf  +  <7l]=-A°R5d(A',)''  +  rf,  (18) 

q—  l  "Jo 

where  Rdd  =  diag  [£)?,*  =  1, 2,  •••,  n0)  and  A°  =  [a1,a2,---,ano]  represent  the  signal 
correlation  matrix  and  the  mixing  matrix  formulated  by  considering  n0  signals  out  of  the 
total  number  of  n  signal  arrivals,  respectively. 

It  is  clear  from  (18)  that  the  SNR  improvement  G  —  L/n0  (we  assume  L  >  n0)  is 
inversely  proportional  to  the  number  of  sources  contributing  to  matrix  D.  Therefore,  from 
the  SNR  perspective,  it  is  best  to  set  n0  =  1,  i.e.,  to  select  the  sets  of  N’  t-f  points  that 
belong  to  individual  signals  one  set  at  a  time,  and  then  separately  evaluate  the  respective 
STFD  matrices. 

This  procedure  is  made  possible  by  the  fact  that  STFD-based  array  processing  is,  in 
essence,  a  discriminatory  technique  in  the  sense  that  it  does  not  require  simultaneous 
localization  and  extraction  of  all  unknown  signals  received  by  the  array.  With  STFDs, 
array  processing  can  be  performed  using  STFDs  of  a  subclass  of  the  impinging  signals  with 
specific  t-f  signatures.  In  this  respect,  the  t-f  based  blind  source  separation  and  direction 
finding  techniques  have  implicit  spatial  filtering,  removing  the  undesired  signals  from  con¬ 
sideration.  It  is  also  important  to  note  that  with  the  ability  to  construct  the  STFD  matrix 
from  one  or  few  signal  arrivals,  the  well  known  m  >  n  condition  on  source  localization 
using  arrays  can  be  relaxed  to  m  >  n„,  i.e.,  we  can  perform  direction  finding  or  source 
separation  with  the  number  of  array  sensors  smaller  than  the  number  of  impinging  signals. 
Further,  from  the  angular  resolution  perspective,  closely  spaced  sources  with  different  t-f 
signatures  can  be  resolved  by  constructing  two  separate  STFDs,  each  corresponding  to  one 
source,  and  then  proceed  with  subspace  decomposition  for  each  STFD  matrix,  followed  by 
an  appropriate  source  localization  method  (MUSIC,  for  example).  The  drawback  using 
different  STFD  matrices  separately  is  of  course  the  need  for  repeated  computations. 
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C.  Signal  and  Noise  Subspaces  Using  STFDs 


The  following  Lemma  provides  the  relationship  between  the  eigendecompositions  of  the 
STFD  matrices  and  the  data  covariance  matrices  used  in  conventional  array  processing. 

Lemma  2  [14]:  Let  \[  >  \°2  >  •  •  •  >  X°no  >  A"o+1  =  A"o+2  =  •■■  =  \°m  =  a  denote 
the  eigenvalues  of  Rxx  =  A°Rdd(A°)H  +  a  I,  which  is  defined  from  a  data  record  of  a 
mixture  of  the  n0  selected  FM  signals.  Denote  the  unit-norm  eigenvectors  associated  with 
A", A"o  by  the  columns  of  S°  =  [s°,...,s®J  ,  and  those  corresponding  to  A°o+1, ...,  \°n 
by  the  columns  of  G®  =  [g°,  ...,g^_nJ.  We  also  denote  A4/  >  A^  >  •  •  •  >  Aj£  >  A ^{+1  = 
^{+2  =  ' ' '  =  >&  =  ”1'  as  the  eigenvalues  of  D  defined  in  (18).  The  superscript  ^  denotes 
that  the  associated  term  is  derived  from  the  STFD  matrix  D.  The  unit-norm  eigenvectors 
associated  with  A4/, ...,  A *Jo  are  represented  by  the  columns  of  =  [s4/, ...,  slJo]  ,  and  those 
corresponding  to  A^{+1,...,A^  are  represented  by  the  columns  of  =  [g4/, ..., gm-nj- 
Then, 

a)  The  signal  and  noise  subspaces  of  and  Gt}  are  the  same  as  S°  and  G°,  respectively. 

b)  The  eigenvalues  have  the  following  relationship: 


( 


(yU  —  a 


i  <  nQ 
n0  <  i  <  m. 


(19) 


An  important  conclusion  from  Lemma  2  is  that,  the  largest  n0  eigenvalues  are  amplified 
using  STFD  analysis.  The  amplification  of  the  largest  nc  eigenvalues  improves  detection 
of  the  number  of  the  impinging  signals  on  the  array,  as  it  widens  the  separation  between 
dominant  and  noise-level  eigenvalues.  Determination  of  the  number  of  signals  is  key  to 
establishing  the  proper  signal  and  noise  subspaces,  and  subsequently  plays  a  fundamental 
role  in  subspace-based  applications.  When  the  input  SNR  is  low,  or  the  signals  are  closely 
spaced,  the  number  of  signals  may  often  be  underdetermined.  When  the  STFD  is  applied, 
the  SNR  threshold  level  and/or  angle  separation  necessary  for  the  correct  determination 
of  the  number  of  signals  are  greatly  reduced. 

Next  we  consider  the  signal  and  noise  subspace  estimates  from  a  finite  number  of  data 
samples.  We  form  the  STFD  matrix  based  on  the  true  ( t ,  /)  points  along  the  IF  of  the  na 
FM  signals. 

Lemma  3  [14],  [10]:  If  the  third-order  derivative  of  the  phase  of  the  FM  signals  is 
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negligible  over  the  time-period  [t  —  L  +  l,t  +  L  —  1],  then 

a)  The  estimation  errors  in  the  signal  vectors  are  asymptotically  (for  N  L)  jointly 
Gaussian  distributed  with  zero  means  and  covariance  matrices  given  by 
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(2i) 

b)  The  orthogonal  projections  of  {g^}  onto  the  column  space  of  are  asymptotically 
(for  N  L)  jointly  Gaussian  distributed  with  zero  means  and  covariance  matrices  given 


E  (stf  (S«)H g^)  (stf  (stf)Hg]fy 

~n0N'\^l{a-X[fySk  ]*< 

=  —  V  —  ~  g)  +  j£<T  0  ,  o\H  r 
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E  (stf  (S*)V)  (s*  (Stf)H  gf  ^  =  0  for  all  *,  j .  (23) 

Prom  (20)-(23),  two  important  observations  are  in  order.  First,  if  the  signals  are  both 
localizable  and  separable  in  the  t-f  domain,  then  the  reduction  of  the  number  of  signals 


59 


from  n  to  n0  greatly  reduces  the  estimation  error,  specifically  when  the  signals  are  closely 
spaced.  The  second  observation  relates  to  SNR  enhancements.  The  above  equations  show 
that  error  reductions  using  STFDs  are  more  pronounced  for  the  cases  of  low  SNR  and/or 
closely  spaced  signals.  It  is  clear  from  (20)-(23)  that,  when  \°k^>  a  for  all  A;  =  1,2, ...,  n0, 
the  results  are  almost  independent  of  L  (suppose  Ar  >■  L  so  that  N'  =  N  —  L+l  ~  N),  and 
therefore  there  would  be  no  obvious  improvement  in  using  the  STFD  over  conventional 
array  processing.  On  the  other  hand,  when  some  of  the  eigenvalues  are  close  to  a  (A£  ~  a, 
for  some  k  =  1, 2, ...,  na),  which  is  the  case  of  weak  or  closely  spaced  signals,  all  the  results 
of  above  three  equations  are  reduced  by  a  factor  of  up  to  G  =  L/n0 ,  respectively.  This 
factor  represents,  in  essence,  the  gain  achieved  from  using  STFD  processing. 

IV.  Blind  Source  Separation 

A.  Source  Separation  Based  on  STFDs 

Blind  source  separation  based  on  STFD  was  first  considered  by  Belouchrani  and  Amin 
[8].  The  first  step  of  STFD-based  blind  source  separation  is  the  whitening  of  the  signal 
part  y (f)  of  the  observation.  This  is  achieved  by  applying  a  whitening  matrix  W  to  y(t), 
i.e.,  an  n  x  N  matrix  satisfying: 

1  T 

lim  -  W.y{t)yH{t)WH  =  WRyyWfi  =  WAAHAH  =  I.  (1) 

T^°°  *  t= i 

WA  is  an  n  x  n  unitary  matrix  U,  and  matrix  A  can  be  written  as 

A  =  W#U  (2) 

where  superscript  *  denotes  pseudo-inverse.  The  whitened  process  z (t)  =  Wx(£)  still 
obeys  a  linear  model, 

z  (t)  =  Wx(t)  =  W  [As  (t)  +  n(t)]  =  Us(£)  +  Wn(i).  (3) 

By  pre-  and  post-multiplying  the  STFD  matrices  Dxx(£, /)  by  W,  we  obtain 

Dzz(£,/)=WDxx(£,/)W".  (4) 

which  is,  in  essence,  the  STFD  of  the  whitened  data  vector  z (£).  From  the  definitions  of 
W  and  U, 

Dzz(t,  /)  =  UDSS(£,  f)UH.  (5) 


60 


Equation  (5)  shows  that  if  Dss(t,  /)  is  diagonal,  which  is  the  case  of  autoterm  points,  then 
any  whitened  data  STFD  matrix  is  diagonal  in  the  basis  of  the  columns  of  the  matrix  U, 
and  the  eigenvalues  of  D zz(t,f)  are  the  diagonal  entries  of  Dss  (£,/).  An  estimate  U  of 
the  unitary  matrix  U  may  be  obtained  as  a  unitary  diagonalizing  matrix  of  a  whitening 
STFD  matrix  for  some  t-f  points  corresponding  to  the  signal  autoterm.  The  source  signals 
can  then  be  estimated  as  s  (t)  —  UWx(t),  and  the  mixing  matrix  A  is  estimated  by 
A  =  W#U. 

In  order  to  reduce  the  noise  effect  as  well  as  the  possibility  of  having  degenerate  eigen¬ 
values  and  subsequently  non-unique  solutions,  the  JD  and  t-f  averaging,  both  discussed  in 
II-C,  can  be  used  to  incorporate  multiple  t-f  points. 

The  method  discussed  above  uses  STFD  matrices  to  estimate  the  unitary  matrix  U, 
but  the  covariance  matrix  is  still  used  for  whitening.  Therefore,  the  advantages  of  STFD 
matrices  are  not  fully  utilized.  Using  the  STFD  matrix  D,  instead  of  the  covariance 
matrix  Rxx,  to  perform  whitening  is  a  reasonable  alternative  [11],  To  avoid  degenerate 
eigenvalues,  the  STFD  matrices  used  for  pre-whitening  and  unitary  matrix  estimation 
should  be  different. 

A. l  Example 

Fig.  1  shows  an  example  of  the  application  of  STFDs  to  the  BSS  problem.  A  three- 
element  equi-spaced  linear  array  is  considered  where  the  interelement  spacing  is  half  a 
wavelength.  Two  chirp  signals  arrive  at  —10°  and  10°,  respectively.  The  number  of  data 
samples  used  to  compute  the  STFD  is  128.  The  number  of  t-f  points  employed  in  the  JD  is 
p=128,  with  equal  number  of  points  on  each  signature.  Fig.  1(b)  shows  the  Choi-Williams 
distributions,  in  which  an  exponential  kernel  is  applied  [16],  of  two  linear  mixtures  of  the 
original  chirp  signals  depicted  in  Fig.  1(a),  corresponding  to  the  data  at  the  first  and  the 
second  sensors.  Using  JD  of  the  STFDs,  we  are  able  to  recover  the  original  signals  from 
their  observed  mixture,  as  shown  in  Fig.  1(c). 

B.  Source  Separation  Based  on  Spatial  Averaging 

Source  separation  based  on  spatial  averaging  is  proposed  by  Mu,  Amin,  and  Zhang  [17]. 
This  method  first  performs  array  averaging  of  the  TFDs  of  the  data  across  the  array, 
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permitting  the  spatial  signature  of  sources  to  play  a  fundamental  role  in  improving  the 
synthesis  performance. 

This  method  is  philosophically  different  from  the  one  described  in  the  previous  section, 
as  it  applies  the  opposite  order  of  operations.  It  first  synthesizes  the  source  signals  from 
the  t-f  domain,  then  proceeds  to  estimate  their  respective  mixing  matrix. 

The  WVD-based  synthesis  techniques  could  be  found  in  [18],  [19].  Herein,  we  apply  the 
method  of  extended  discrete-time  Wigner  distribution  (EDTWD),  introduced  in  [19],  to 
the  output  of  array  averaged  WVD.  The  advantage  of  using  the  EDTWD  lies  in  the  fact 
that  it  does  not  require  a  priori  knowledge  of  the  source  waveform,  and  thereby  avoids  the 
problem  of  matching  the  two  “uncoupled”  vectors  (even-indexed  and  odd-indexed  vectors). 
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The  overall  synthesis  procedure  is  summarized  in  the  following  steps. 

1.  Given  the  received  data  of  the  ith  sensor  xt(t),  compute  the  EDTWD 

W XiXi(t,  f)=  £  Xi(t  +  k/2)x*{t  -  k/2)e~^, 

k:{t+k/2)eZ 

t  =  0,±0.5,±1,....  (6) 


2.  Apply  the  averaging  process,  that  is,  summing  the  EDTWD  across  the  array 

i  m 

W(t,f)  =  (7) 

m  k= 1 

3.  Place  an  appropriate  t-f  mask  on  W (t,  /)  such  that  only  the  desired  signal  autoterms 
are  retained. 

4.  Take  the  IFFT  of  the  masked  WVD  W(t,  /) 


1 

5.  Construct  the  matrix  Q  =  [qu\  with 


(8) 

(9) 


6.  Apply  eigendecomposition  to  the  Hermitian  matrix  [Q  +  Q^]  and  obtain  the  maximum 
eigenvalue  Amax  and  the  associated  eigenvector  u.  The  desired  signal  is  given  by 


S  opt  =  eja 


(10) 


where  a  is  an  unknown  value  representing  the  phase. 

7.  Repeat  step  3  through  6  until  all  source  signals  di(t),d,2(t), . . . ,  d^t)  are  retrieved. 

The  averaging  in  step  2  mitigates  the  crossterms  and  enforces  the  autoterms.  As  such, 
the  source  t-f  signatures  become  easier  to  identify,  mask,  and  synthesize.  It  is  notewor¬ 
thy  that  (7)  will  completely  suppress  the  cross-TFDs  for  sources  with  orthogonal  spatial 
signatures. 

Upon  synthesizing  all  the  source  signals,  we  could  utilize  these  signal  waveforms  to 
estimate  the  mixing,  or  array,  matrix  A  through  the  minimization  of  the  mean  square 
error  (MSE), 

e  =  X  llx(0  -  A§(<)ll2  -  (11) 

t= l 
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This  results  in 


A  =  rRT1.  (12) 

where  R  =  J2t= 1  d(t)dH(t)  represents  the  estimated  signal  source  covariance  matrix,  and 
r  =  with  fj  =  J2t=i  x(f)d*  (i)  being  the  correlation  vector  between  the  data 

vector  received  across  the  array  and  the  tth  source  signal  di(t). 

B.l  Example 

We  consider  three  parallel  chirp  signals.  The  signals  arrive  with  DOAs  of  —20°,  0°  and 
20°,  with  the  respective  start  and  end  frequencies  given  by  (0.97T,  0.57t),  (0.667T,  0.267t), 
and  (0.57T,  0.l7r),  respectively.  The  length  of  the  signal  sequence  is  set  to  N  =  128.  The 
input  SNR  is  —5  dB.  The  crossterm  of  di(t)  and  d3(t)  lies  closely  to  the  t-f  signature  of 
d2{t). 

We  first  consider  the  single  antenna  case.  Fig.  2  depicts  both  the  WVD  of  the  signal 
arrival  and  the  WVD  of  the  synthesized  d2(t).  The  signal  is  significantly  corrupted  by  the 
crossterm  of  di(t)  and  d3(t)  as  well  as  by  the  noise  components. 

Figure  3  depicts,  for  16  sensor  scenarios,  the  array  averaged  WVD  and  the  respective 
d2(t).  Upon  averaging,  both  noise  and  crossterms  are  sufficiently  reduced  to  clearly  man¬ 
ifest  the  individual  source  t-f  signatures.  The  signals  could,  therefore,  be  individually 
recovered  by  placing  appropriate  masks  in  the  t-f  region.  The  significance  of  using  array 
sensors  is  evident  in  Fig.  3. 


Averaged  WVD  from  1  sensors  WVD  of  synthesized  signal 


time  time 

Fig.  2.  WVD  and  the  synthesized  signal  (M  =  1). 
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WVD  of  synthesized  signal 


Averaged  WVD  from  16  sensors 


Fig.  3.  WVD  and  the  synthesized  signal  (M  =  16). 


C.  Effect  of  Crossterms  between  Source  Signals 

In  this  section,  we  examine  the  effect  of  the  t-f  crossterms  on  source  separation  perfor¬ 
mance  [20].  To  simplify  the  problem,  we  assume  that  Rdd  is  an  identity  matrix.  When 
crossterms  are  present  at  the  off-diagonal  elements  of  the  TFD  matrix  Ddd(t,  /),  then 

Ddd  (t,  f )  =  P(f,  f)G(t,  f)PH(t,  f )  (13) 

where  G (i,  /)  is  the  diagonal  matrix  with  the  eigenvalues  at  the  diagonal  elements,  and 
P(f,  /)  is  the  matrix  whose  columns  are  the  corresponding  eigenvectors.  Note  that  all  the 
above  matrices  depend  on  the  selected  (t,  /)  point.  From  (13),  the  STFD  matrix  of  the 
data  vector  under  noise-free  conditions  becomes 

Dxx(f,  /)  -  ADss(f,  /) Ah  =  AP(t,  f)G(t,  f)PH(t ,  f)AH  (14) 

and  the  STFD  matrix  of  the  whitened  array  signal  vector  is 

DaI(f ,  /)  =  WAP(i,  f)G(t,  f)PH(t,  1)  A*W“  (15) 

Since  G(t,  /)  is  diagonal,  WAP(i,  /)  is  unitary.  Therefore,  the  source  separation  method 
will  assume  WAP(t,  /)  as  the  unitary  matrix  and  estimates  the  mixing  matrix  as 

A  =  W#WAP(t,  /)  =  AP(f,  /)  (16) 

which  is  dependent  on  the  unitary  matrix  P (£,/).  Furthermore, 

A*  A  =  [AP(<, /)}*  A  =  P  «(t, !).  (17) 
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Matrix  A  should  be  close  to  the  true  one  A  so  that  A* A  well  approximates  the  iden¬ 
tity  matrix.  The  following  variable  measures,  at  separations,  the  ratio  of  the  power  of 
interference  of  the  qth.  source  signal  to  the  power  of  the  pth  source  signal  [8] 

/m  =  e|(a#a)J\  (18) 

where  (A#A)  denotes  the  (p,q) th  element  of  matrix  A* A.  The  following  global 
rejection  level  is  used  to  evaluate  the  overall  performance  of  a  blind  source  separation 
system 

Iperf  =  52  IpQ'  (19) 

qjp 

For  the  mixing  matrix  estimation  given  in  (17),  the  global  rejection  level  is  approximated 
by  the  following  normalized  global  rejection  level  [20] 

t'perf  =  [diagonal  (A#A)]  A* A  =  Y  IPto(*>  /)l~2  ~  n •  (20) 

9=1 

where  diagonal( F)  denotes  the  matrix  formed  by  the  diagonal  elements  of  F.  In  general, 
since  the  absolute  values  of  pqq(t,  f)  are  always  equal  to  or  smaller  than  1,  the  global 
rejection  level  Iperf  takes  a  positive  value.  It  is  clear  that  Iperf  =  0  only  when  pqq(t ,  /)  =  1 
holds  true  for  all  q.  That  is,  P  is  an  identity  matrix,  which  implies  that  there  is  no 
off-diagonal  non-zero  elements  in  matrix  D<id(i,  /),  i.e.,  no  crossterms. 

D.  Source  Separation  Based  on  Joint  Diagonalization  and  Joint  Anti-Diagonalization 

From  the  previous  subsection,  it  is  clear  that  care  must  be  exercised  when  dealing  with 
crossterms.  The  method,  proposed  by  Belouchrani,  Abed-Meraim,  Amin,  and  Zoubir  [21], 
carefully  and  properly  exploits  both  autoterms  and  crossterms  of  the  TFDs  for  improved 
source  separation.  This  approach  is  based  on  the  simultaneous  diagonalization  and  anti- 
diagonalization  of  a  combined  set  of  autoterm  and  crossterm  TFD  matrices,  respectively. 
The  auto-STFD  and  cross-STFD  are  defined  as 

D“s(i,  /)  =  Dss(t,  /)  for  autoterm  t-f  points  (21) 

and 

Dcss(t,  f)  =  Dss(t,  /)  for  crossterm  t-f  points.  (22) 
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Since  the  off-diagonal  elements  of  Dss(t,  /)  are  crossterms,  the  auto-STFD  matrix  is  quasi 
diagonal  for  each  t-f  point  that  corresponds  to  a  true  power  concentration,  i.e.  signal 
autoterm.  Similarly,  since  the  diagonal  elements  of  Dss(t,/)  are  auto-terms,  the  cross- 
STFD  matrix  is  quasi  anti-diagonal  (i.e.  its  diagonal  entries  are  close  to  zero)  for  each  t-f 
point  that  corresponds  to  a  crossterm. 

As  discussed  earlier,  JD  can  be  used  to  incorporate  multiple  autoterm  t-f  points.  Simi¬ 
larly,  the  joint  anti-diagonalization  ( JAD)  is  appropriate  to  incorporate  multiple  crossterm 
t-f  points.  By  selecting  crossterm  t-f  points,  the  data  cross-STFD  will  have  the  following 
structure, 

!£,(*,/)  =  UDJ,(t,/)UH  (23) 

where  D £s(t,/)  is  anti-diagonal.  The  JAD  searches  for  the  unitary  matrix  that  anti¬ 
diagonalizes  a  combined  set  {D £x(t;,  fi)\i  =  1,  •  •  • ,  q)  of  q  STFD  matrices.  The  procedure 
for  anti-diagonalization  of  a  single  m  x  m  matrix  N  is  explained  in  [21]  and  is  equivalent 
to  the  maximization  of  the  criterion 

m  ~ 

<7(N,V)^f-£|vfNvi|  (24) 

i=l 

over  the  set  of  unitary  matrices  V  =  [vi,  •  •  • ,  vm]. 

The  combined  JD  and  JAD  of  two  sets  =  l..p}  and  {Nk\  k  =  l..q}  ofmxm 

matrices  is  defined  as  the  maximization  of  the  JD/JAD  criterion: 

C(V)  =  e(e  IvfM.vJ2  -  ±  |v"Ntvi|2)  (25) 

1=1  \k= 1  Jfc— 1  / 

over  the  set  of  unitary  matrices  V  =  [vl5  •  •  • ,  vm]. 

D.l  Selection  of  Autoterm  and  Crossterm  Points 

The  success  of  the  JD  or  JAD  of  STFD  matrices  in  determining  the  unitary  matrix  U 
depends  strongly  on  the  correct  selection  of  the  autoterm  and  crossterm  points.  Therefore, 
it  is  crucial  to  have  a  selection  procedure  that  is  able  to  distinguish  between  autoterm  and 
crossterm  points  based  only  on  the  STFD  matrices  of  the  observation.  A  selection  approach 
was  proposed  in  [21]  to  exploit  the  anti-diagonal  structure  of  the  crossterm  STFD  matrices. 
More  precisely, 

Trace(D^x(t,  /))  =  Trace(UD‘s(t,  f)VH)  =  TVace(Dscs(t,  /))  «  0. 


67 


Based  on  this  observation,  the  following  testing  procedure  can  be  defined: 

if  Trace (Dxx (t ,  f  ) )  <  f  — »  decide  that  (£,  /)  is  a  crossterm 
norm(Dxx(t,  /)) 

if  Trace(Dxx(t,  /))  — *  decide  that  (t,  /)  is  an  autoterm 

norm(Dxx  (£,/))  V 

where  e  is  a  ‘small’  positive  real  scalar. 

D.2  Example 

We  consider  a  uniform  linear  array  of  m  =  3  sensors  having  half  wavelength  spacing 
and  receiving  signals  from  n  =  2  sources  in  the  presence  of  white  Gaussian  noise.  The 
sources  arrive  from  different  directions  0\  =  10  and  82  —  20  degrees.  The  emitted  signals 
are  two  chirps.  The  WVD  is  computed  over  1024  samples  and  eight  STFD  matrices  are 
considered. 

We  compare  in  Fig.  4  the  performance  of  the  JD-based  algorithm,  introduced  in  IV-A, 
and  the  JD/JAD  algorithm,  for  SNRs  in  the  range  [5  -  20  dB].  The  mean  rejection  levels 
are  evaluated  over  100  Monte  Carlo  runs.  In  this  case,  the  new  algorithm  performs  slightly 
better  than  the  JD-based  algorithm. 


-30-  'V 

-31  - 

-32' - - - » - 1 - 

5  10  IS  20 

SNR  (dB) 

Fig.  4.  Mean  rejection  level  versus  input  SNR  for  JD  and  JD/AJD  based  source  separation  methods. 

V.  Direction  Finding 

A.  Time-Frequency  MUSIC 

The  t-f  MUSIC  was  proposed  by  Belouchrani  and  Amin  [9],  and  its  performance  is 
analyzed  by  Zhang,  Mu,  and  Amin  [14].  Without  loss  of  generality,  we  consider  one- 
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dimensional  direction  finding  where  the  DOAs  are  described  by  9.  First,  recall  that  the 
DOAs  are  estimated  in  the  MUSIC  technique  by  determining  the  n  values  of  9  for  which 
the  following  spatial  spectrum  is  maximized  [22], 


JW(«)  =  [>'(#>)] =  [a* (0)  (i  -  SS")  a(0)] (1) 

where  a(9)  is  the  steering  vector  corresponds  to  9.  The  variance  of  those  estimates  in  the 
MUSIC  technique,  assuming  white  noise  processes,  is  given  by  [15] 


E  (u h  -  u >i) 


2 _  1  aff(^)Ua(^) 

2  N  h(Qi) 


where  a;,-  =  (27rd/A)  sin  9%  is  the  spatial  frequency  associated  with  DOA  9U  and  a),  is  its 
estimate  obtained  from  the  MUSIC.  Moreover,  U  is  defined  in  (11),  and 

h(9i)  =  dH(9i)GGHd(9i),  with  d(0,)  =  da(9i)/du>.  (3) 

Similarly,  for  t-f  MUSIC  with  nQ  signals  selected,  the  DOAs  are  determined  by  locating 
the  nQ  peaks  of  the  spatial  spectrum  defined  from  the  n0  signals’  t-f  regions, 

filvV)  =  [a"(#)e"  (G")"  a(«)] =  [a"(»)  (i  -  S1^  (S1')*)  a(0)]  .  (4) 

GU  and  can  be  obtained  by  using  either  JBD  or  t-f  averaging  (section  II-C).  When 
the  t-f  averaging  is  used,  using  the  results  of  Lemmas  2  and  3,  the  variance  of  the  DOA 
estimates  based  on  t-f  MUSIC  is  obtained  as  [14], 

E  (j\tf  .  1  aif(«i)U‘/a(«i) 

-w0  -w  vm  (5) 


where  is  the  estimate  of  a jt,  is  defined  in  (22),  and 

htf(9i)  =  dH(9t)Gtf  (Gtf)Hd &). 
Note  that  htf  {9)  =  h(9i)  if  n0  —  n. 


A.l  Examples 

Consider  a  uniform  linear  array  of  8  sensors  spaced  by  half  a  wavelength,  and  an  obser¬ 
vation  period  of  1024  samples.  Two  chirp  signals  emitted  from  two  sources  positioned  at 


angle  and  02-  The  start  and  end  frequencies  of  the  signal  source  at  61  are  usi  =  0  and 
toe\  =  7r,  while  the  corresponding  two  frequencies  for  the  other  source  at  02  are  los2  =  ^ r 
and  u!e2  =  0,  respectively. 

Fig.  5  displays  the  variance  of  the  estimated  DOA  6\  versus  SNR  for  the  case  (#1,  $2)  = 
(—10°,  10°).  The  curves  in  this  figure  show  the  theoretical  and  experimental  results  of  the 
conventional  MUSIC  and  t-f  MUSIC  (for  L= 33  and  129).  The  Cramer-Rao  bound  (CRB) 
is  also  shown  in  Fig.  5  for  comparison.  Both  signals  were  selected  when  performing  t-f 
MUSIC  (n0  =  n  =  2).  Simulation  results  were  averaged  over  100  independent  Monte- 
Carlo  runs.  The  advantages  of  t-f  MUSIC  in  low  SNR  cases  are  evident  from  this  figure. 
The  experiment  results  deviate  from  the  theoretical  results  for  low  SNR,  since  we  only 
considered  the  lowest  order  of  the  coefficients  of  the  perturbation  expansion  in  deriving 
the  theoretical  results  [14].  Fig.  6  shows  estimated  spatial  spectra  at  SNR=— 20  dB  based 
on  t-f  MUSIC  ( L  =  129)  and  the  conventional  MUSIC.  The  t-f  MUSIC  spectral  peaks  are 
clearly  resolved. 

Fig.  7  shows  examples  of  the  estimated  spatial  spectrum  based  on  t-f  MUSIC  ( L  —  129) 
and  the  conventional  MUSIC  where  the  angle  separation  is  small  (Q\  =  —2.5°,  d2  —  2.5°). 
The  input  SNR  is  —5  dB.  Two  t-f  MUSIC  algorithms  are  performed  using  two  sets  of 
t-f  points,  each  set  belongs  to  the  t-f  signature  of  one  source  (n0  =  1).  It  is  evident 
that  the  two  signals  cannot  be  resolved  when  the  conventional  MUSIC  is  applied,  whereas 
by  utilizing  the  signals’  distinct  t-f  signatures  and  applying  t-f  MUSIC  separately  for 
each  signal,  the  two  signals  become  clearly  separated  and  reasonable  DOA  estimation  is 
achieved.  It  is  noted  that  there  is  a  small  bias  in  the  estimates  of  t-f  MUSIC  due  to  the 
imperfect  separation  of  the  two  signals  in  the  t-f  domain. 

B.  Time- Frequency  Maximum  Likelihood  Method 

In  this  section,  we  introduce  the  time-frequency  maximum  likelihood  (t-f  ML)  methods. 
This  method  was  proposed  by  Zhang,  Mu,  and  Amin  [10]  to  deal  with  coherent  nonsta¬ 
tionary  sources.  For  conventional  ML  methods,  the  joint  density  function  of  the  sampled 
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SNR  <d8) 

Fig.  5.  Variance  of  DOA  estimation  versus  input  SNR. 


t-f  MUSIC 
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Fig.  6.  Estimated  spatial  spectra  of  MUSIC  and  t-f  MUSIC. 


data  vectors  x(l),x(2),  ...,x(Ar),  is  given  by  [23] 

/(x(l),...,x(A)) 

=  n  ^p]exP  (~a  [x(i)  -  Ad(i)]H  [x(z)  -  Ad(i)0  ’ 

where  det[-]  denotes  the  matrix  determinant.  It  follows  from  (7)  that  the  log-likelihood 
function  of  the  observations  x(l),x(2),  . . .  ,  x(/V),  is  given  by 

1  N 

L  =  —mN\na - [x(i)  —  Ad(i)]H  [x(i)  —  Ad(i)] .  (8) 


To  carry  out  this  minimization,  we  fix  A  and  minimize  (8)  with  respect  to  d.  This  yields 


the  well-known  solution 


d(i)  =  [A*  A]  1  A  Hx(i). 
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t-f  MUSIC 


Fig.  7.  Estimated  spatial  spectra  of  MUSIC  and  t-f  MUSIC  for  closely  spaced  signals. 

We  can  obtain  the  concentrated  likelihood  function  as  [23] 

FMl(Q)  =  tr  { [I  -  A(Ah A)-1  Ah]  Rxx}  ,  (10) 

where  tr( A)  denotes  the  trace  of  A.  The  ML  estimate  of  ©  is  obtained  as  the  minimizer 
of  (10).  Let  u>i  and  u>i,  respectively,  denote  the  spatial  frequency  and  its  ML  estimate 
associated  with  0i,  then  the  estimation  error  (cot  —  Lot)  are  asymptotically  (for  large  N ) 
jointly  Gaussian  distributed  with  zero  means  and  the  covariance  matrix  [24] 

E  [(A  -  a,,)2]  =  ^  [Re(H  ©  Rjd)] 

x  Re  [H  ©  (Rdd  A"UARdd)1']  [Re(H  ©  Rjd)] 
where  ©  denotes  Hadamard  product,  U  is  defined  in  (11).  Moreover, 

H  =  CH  [i  -  A(Ai/A)“1Aff]  C,  with  C  =  dA/du>.  (12) 

Next  we  consider  the  t-f  ML  method.  As  we  discussed  in  the  previous  section,  we  select 
n0  <  n  signals  in  the  t-f  domain.  The  concentrated  likelihood  function  defined  from  the 
STFD  matrix  is  similar  to  (10)  and  is  obtained  by  replacing  Rxx  by  D, 

=  tr  [i  -  A0  ((A°)"A°)-1  (A*)"]  D.  (13) 

Therefore,  the  estimation  error  (T  ~  T  associated  with  the  t-f  ML  method  are  asymp¬ 
totically  (for  N  L)  jointly  Gaussian  distributed  with  zero  means  and  the  covariance 
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matrix  [10] 


E  [(<&!'  -  a,,)2]  =  ^  [Re(H"  0  Ddd)]_1 
xRe  [h”0  (Ddd(A”)"U"A“Dd<1)T]  [Re(H“  0  Djd)]"‘ 

=  2^[fe(H”0(R2-)7T' 

xRe  [h“  ©  (Rdd(A°)ffU^A°Rdd)T]  [Re  ((H°  ©  RSd)T)] , 
where  is  defined  in  (22),  and 

H°  =  (C°)ff  [i- A°((A0)ffAo)_1(Ao)H|  C°,  with  C°  =  dA°/du.  (15) 

In  the  case  of  n0  =  n ,  then  H°  =  H,  and  C°  =  C. 

The  signal  localization  in  the  t-f  domain  enables  us  to  select  fewer  signal  arrivals.  This 
fact  is  not  only  important  in  improving  the  estimation  performance,  particularly  when 
the  signals  are  closely  spaced,  but  also  reduces  the  dimension  of  the  optimization  problem 
solved  by  the  maximum  likelihood  algorithm,  and  subsequently  reduces  the  computational 
requirement. 

B.l  Examples 

To  demonstrate  the  advantages  of  t-f  ML  over  both  the  conventional  ML  and  the  t- 
f  MUSIC,  consider  a  uniform  linear  array  of  8  sensors  separated  by  half  a  wavelength. 
Two  FM  signals  arrive  from  (0i,d2)  =  (—10°,  10°)  with  the  IFs  fi(t )  =  0.2  +  O.lt/N  + 
0.2sin(27rt/AO  and  f2(t)  =  0.2  +  0.1t/N  +  0.2sin(27rt/iV  +  tt/2),*  =  1, ...,  N.  The  SNR  of 
both  signals  is  —20  dB,  and  the  number  of  snapshots  used  in  the  simulation  is  N  =  1024. 
We  used  L= 129  for  t-f  ML.  Fig.  8  shows  (8\,02)  that  yield  the  minimum  values  of  the 
likelihood  function  of  the  t-f  ML  and  the  ML  methods  for  20  independent  trials.  It  is 
evident  that  the  t-f  ML  provides  much  improved  DOA  estimation  over  the  conventional 
ML. 

In  the  next  example,  the  t-f  ML  and  the  t-f  MUSIC  are  compared  for  coherent  sources. 
The  two  coherent  FM  signals  have  common  IFs  —  0.2+0.U/N+0.2 sin(2nt/N),  t  = 

1, ...,  iV,  with  7r/2  phase  difference.  The  signals  arrive  at  (61, 02)  =  (—2°,  2°).  The  SNR  of 
both  signals  is  5  dB  and  the  number  of  snapshots  is  1024.  Fig.  9  shows  the  contour  plots 
of  the  likelihood  function  of  the  t-f  ML  and  the  estimated  spectra  of  t-f  MUSIC  for  three 
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independent  trials.  It  is  clear  that  the  t-f  ML  can  separate  the  two  signals  whereas  the  t-f 
MUSIC  cannot. 


1-,. 


Fig.  8.  {0\,  Of)  which  minimize  the  t-f  ML  (upper)  and  ML  (lower)  likelihood  functions. 

C.  Effect  of  Crossterms 

Identifying  the  sources’  t-f  signatures  often  require  searching  the  t-f  domain  for  peak 
values.  In. some  cases,  these  values  correspond  to  crossterms.  Building  the  STFDs  around 
only  crossterms  or  a  mixture  of  autoterm  and  crossterms  and  its  effect  on  the  t-f  MUSIC 
performance  is  considered  by  Amin  and  Zhang  [25].  To  understand  the  tole  of  crossterms 
in  direction  finding,  it  is  important  to  compare  the  crossterms  to  the  cross-correlation 
between  signals  in  conventional  array  processing,  whose  properties  are  familiar.  The  source 
TFD  matrix  takes  the  following  general  form, 

' Ddldl(tJ )  Ddld2(tJ)  •••  Ddldn(tJ) 

Dd2dl(t,  f)  Dd2d2(t,  f)  *••  Dd2dn(t,  /) 

°dd  (tj)=  .  .  ...  . 

- Ddndl(tJ )  Ddnd2(tJ)  •••  Ddndn{t,f) 

On  the  other  hand,  the  covariance  matrix  of  correlated  source  signals  is  given  at  the 
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MML 


Fig.  9.  Contour  plots  of  t-f  ML  likelihood  function  (upper)  and  spatial  spectra  of  t-f  MUSIC  (lower). 


form 


R 


d\ d\ 


Rdd  = 


Rd\d2  *  *  *  Rd\dn 
Rd2d2  *  *  *  Rd2dn 


\~Rdndi  Rdnd2 


R 


'dndn  J 


(17) 


where  the  off-diagonal  element  R^dj  =  E[di(t)d*{t)]  represents  the  correlation  between 
source  signals  di  and  dj.  Direction  finding  problems  can  usually  be  solved  when  the  signals 
are  partially  correlated,  however,  full  rank  property  of  the  source  covariance  matrix  Rdd 
is  a  necessary  condition. 

Comparing  equations  (16)  and  (17),  it  is  clear  that  the  cross-correlation  terms  and  the 
crossterms  have  analogous  forms.  However,  the  correlation  matrix  in  (17)  is  defined  for 
stationary  signal  environments,  whereas  the  source  TFD  matrix  in  (16)  is  defined  at  a 
(t,  /)  point  and  its  value  usually  varies  with  respect  to  t  and  /.  Detailed  observations  are 
made  through  the  following  example. 


C.l  Example 

Consider  a  six-element  linear  array  with  half- wavelength  inter-element  spacing,  and  two 
chirp  signal  arrival.  The  start  and  end  frequencies  of  the  first  signal  d\(t)  are  /is  =  0.1 
and  /ie  =  0.5,  and  those  for  the  second  signal  d2(t)  are  }2s  =  0  and  /2e  =  0.4,  respectively. 
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The  SNR  is  10  dB  for  each  signal,  and  the  DOAs  of  the  two  signals  are  61  =  —5°  and 
62  =  5°,  respectively.  The  number  of  samples  is  256.  PWVD  is  used  and  the  window 
length  is  N  =  129. 

We  consider  the  autorerms  and  the  crossterms  over  the  following  two  regions:  i)  au¬ 
toterm  regions  ( f ,  /1)  with  /i(f)  =  0.1  +  OAt/N  and  (f,  /2)  with  /2(f)  =  OAt/N,  where  the 
autoterms  are  dominant;  and  ii)  crossterm  region  (t,  fc)  with  /c(f)  =  [fi(t)  +  /2(f)]/2  — 
0.05  +  OAt/N,  where  the  crossterms  are  dominant.  Both  the  autoterm  and  crossterm 
regions  have  large  peak  values  and  are  most  likely  to  be  selected. 

i)  Autoterm  regions.  In  the  autoterm  region  of  di(t),  (t,  /1),  the  autoterm  of  di(t)  is 
constant.  The  autoterm  of  d2(f)  and  the  crossterm  between  d\(t)  and  d2(t)  are  relatively 
small.  Since  the  STFD  matrix  in  the  autoterm  region  has  dominant  diagonal  elements  with 
constant  values,  incorporating  only  autoterm  points,  either  by  JBD  or  by  t-f  averaging, 
usually  provides  good  direction  finding  performance. 

ii)  Crossterm  regions.  In  this  region  the  crossterms  Ddld2(t,  f)  =  Dd2dl(t,f)  are  domi¬ 
nant.  Therefore,  the  source  TFD  matrix  on  the  crossterm  region  is  nearly  anti-diagonal. 
Note  that  this  source  TFD  matrix  is  still  full  rank  (although  not  necessarily  positive  def¬ 
inite).  Accordingly,  the  noise  subspace  can  be  properly  estimated,  even  when  only  the 
crossterm  points  are  selected.  However,  since  the  crossterms  change  with  time  t ,  taking 
both  positive  and  negative  values,  summing  them  at  different  (f, /)  points  yields  small 
smoothed  values.  Therefore,  the  t-f  averaging  is  expected  to  yield  degraded  performance 
in  some  cases.  Performing  JBD  instead  of  t-f  averaging  avoids  such  risk. 

Table  I  shows  the  DOA  variance  of  signal  di(t)  obtained  from  100  independent  Monte- 
Carlo  runs  of  t-f  MUSIC.  Both  JBD  and  t-f  averaging  are  considered  for  four  cases,  namely, 
(a)  autoterm  regions  /(f)  =  fi(t)  and  f(t)  =  f2(t),  (b)  crossterm  region  f(t )  =  [/i(t)  + 
/2(f)]/2,  (c)  autoterm  and  crossterm  regions  /(f)  =  /i(f),  /(f)  =  /2(f),  and  /(f)  =  [/1(f)  + 
/2(f)]/2,  and  (d)  autoterm  region  of  the  first  signal,  /(f)  =  /i(f).  Although  both  the 
JBD  and  t-f  averaging  resolve  the  signals  in  all  the  four  cases,  it  is  evident  that  the  JBD 
outperforms  the  t-f  averaging,  particularly  when  the  crossterm  region  is  involved.  Case 
(d)  in  which  only  one  of  the  two  signals  is  selected  has  the  best  performance  for  both 
methods  of  JBD  and  t-f  averaging.  An  interesting  observation  is  that,  in  case  (b),  where 
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only  the  crossterm  region  is  used,  the  JBD  yields  second  best  performance,  whereas  the 
t-f  averaging  shows  its  worst  performance. 

TABLE  I 

Variances  of  DOA  estimates 


case  (a) 

case  (b) 

case  (c) 

case  (d) 

JBD 

0.156° 

0.154° 

0.180° 

O 

FO 

A 

0 

T-f  averaging 

0.179° 

0.339° 

0.199° 

0.161° 

VI.  Spatial  Ambiguity  Functions 

The  spatial  ambiguity  function  is  proposed  by  Amin,  Belouchrani,  and  Zhang  [26].  The 
discrete  form  spatial  ambiguity  function  (SAF)  matrix  of  a  signal  vector  x(f)  is  defined  as 

OO 

Dxx(0,t)=  Y,  x(u  +  r/2)xff  (u  -  r/2)e^“  (1) 

U—  —  OO 

where  6  and  r  are  the  frequency-lag  and  the  time-lag,  respectively.  In  noise-free  environ¬ 
ment,  x(t)  =  Ad(t).  In  this  case, 

Dxx(0,r)  =  ADdd(0,T)A*.  (2) 

Equation  (2)  is  similar  to  the  formula  that  has  been  commonly  used  in  blind  source 
separation  and  DOA  estimation  problems,  relating  the  data  correlation  matrix  to  the 
signal  correlation  matrix.  Here,  these  matrices  are  replaced  by  the  data  SAF  and  signal 
ambiguity  function  matrices,  respectively.  The  two  subspaces  spanned  by  the  principle 
eigenvectors  of  Dxx(0,t)  and  the  columns  of  A  are  identical.  This  implies  that  array 
signal  processing  problems  can  be  approached  and  solved  based  on  the  SAF. 

By  replacing  the  STFD  matrix  Dxx(t,  /)  by  the  SAF  matrix  D ^(d,  r),  we  can  easily  de¬ 
rive  the  ambiguity-domain  source  separation  methods  and  the  ambiguity-domain  MUSIC 
(AD  MUSIC)  [26],  following  the  same  procedures  described  in  IV  and  V. 

The  SAFs  have  the  following  two  important  offerings  that  distinguish  them  from  other 
array  spatial  functions. 

1)  The  crossterms  in  between  source  signals  reside  on  the  off-diagonal  entries  of  matrix 
Ddd($,  t),  violating  its  diagonal  structure,  which  is  necessary  to  perform  blind  source 
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separation.  In  the  ambiguity  domain,  the  signal  autoterms  are  positioned  near  and  at  the 
origin,  making  it  easier  to  leave  out  crossterms  from  matrix  construction. 

2)  The  autoterms  of  all  narrowband  signals,  regardless  of  their  frequencies  and  phases, 
fall  on  the  time-lag  axis  (0  =  0),  while  those  of  the  wideband  signals  fall  on  a  different  (0,  r) 
region  or  spread  over  the  entire  ambiguity  domain.  Therefore,  the  SAF  is  a  natural  choice 
for  recovering  and  spatially  localizing  narrowband  sources  in  broadband  signal  platforms. 

VII.  Sensor- Angle  Distributions 

In  this  section  we  use  quadratic  distributions  to  address  the  problem  of  characterizing 
the  power  attributed  to  near-field  scattering  local  to  an  array  of  sensors.  The  proposed 
method  is  based  on  the  quadratic  sensor-angle  distribution  (SAD),  previously  called  the 
spatial  Wigner  distribution  [27].  This  distribution  is  a  characterization  of  the  power  at 
every  angle  for  each  sensor  in  the  array.  It  is  altogether  different  than  the  STFDs  discussed 
in  the  previous  sections.  These  two  types  of  distributions  have  different  structures  and 
objectives. 

In  the  SAD,  near-held  sources  have  different  angles  for  the  various  array  sensors.  The 
SAD  is  a  joint-variable  distribution  and  a  dual  in  sensor  number  and  angle  to  Cohen’s 
class  of  TFDs  [1].  We  use  a  known  test  source  to  illuminate  the  local  scatterer  distribution 
we  wish  to  characterize.  Orthogonal  subspace  projection  techniques  are  then  applied  to 
the  array  data  to  suppress  the  direct  propagation  path  from  the  test  source  so  as  to  reveal 
the  less  powerful  local  scatter.  An  example  from  the  area  of  high-frequency  surface  wave 
radar  is  provided  for  illustration. 

A  typical  surface- wave  radar  receiving  array  may  consist  of  between  8  and  64  sensors  and 
can  be  hundreds  of  meters  or  indeed  more  than  1km  in  total  length.  It  is  typically  sited 
on  a  coastal  beach  which  may  or  may  not  provide  a  uniform  transition  from  land  to  sea. 
The  coast  may  in  fact  be  a  bay  in  which  case  the  land  sea  boundaries  beyond  either  end  of 
the  array  may  cause  near-field  scattering  and  distort  the  wave-front  arriving  at  the  array. 
There  may  be  other  locally  sited  structures,  such  as  buildings  and  fences,  which  can  be 
the  source  of  local  scatter  (consider  that  the  wavelength  of  the  radar  signal  is  between  30- 
100m).  This  makes  achieving  very  low  sidelobe  spatial  beams  with  a  classical  beamformer 
a  difficult  problem  and  can  render  the  receiver  system  vulnerable  to  interference  through 
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beam  sidelobes  (possibly  via  sky  wave  propagation). 

The  near-field  scatter  produced  by  these  mechanisms  are  correlated  with  the  desired 
direct  far-field  radar  return  from  targets  (and  clutter).  This  scatter  is  typically  approx¬ 
imately  20-40dB  weaker  than  the  direct  signal.  Without  compensating  for  the  effects 
of  local  scatter,  it  is  possible  to  achieve  classical  beam  sidelobes  of  30-35dB,  however  in 
general  the  remaining  components  of  the  receiving  system  can  sustain  substantially  higher 
performance  [28]. 

The  effect  of  local  scatter  on  beamforming  must  be  mitigated  in  order  for  the  radar 
system  to  realize  the  inherent  sidelobe  capability  as  set  by  the  radar  equipment  (as  distinct 
from  the  sensors)  [28].  A  first  step  to  achieving  this  is  to  characterize  the  local  scatter 
distribution.  A  means  of  performing  this  characterization  using  techniques  derived  from 
time-frequency  analysis  is  the  focus  of  the  remaining  sections  of  this  chapter. 

A  generalization  of  the  spatial  Wigner  distribution  introduced  in  [27]  is  provided  and 
combined  with  orthogonal  projection  techniques  for  detection,  classification,  and  charac¬ 
terization  of  near-field  and  far-field  sources  lying  in  the  field  of  view  of  the  multi-antenna 
receiver. 

A.  Signal  Model 
A.l  Geometry 

Consider  a  linear  equi-spaced  array  of  M  sensors  placed  on  a  flat  plane  in  a  two  dimen¬ 
sional  surface.  Assume  that  sensor  position  errors  are  negligible  and  the  gain  and  phase 
of  all  sensors  and  corresponding  data  acquisition  equipment  are  accurately  matched.  It 
is  also  assumed  that  the  array  is  narrowband,  i.e.,  the  reciprocal  of  the  bandwidth  of 
any  signals  received  is  large  compared  with  the  propagation  delay  across  the  array.  The 
wavelength  of  all  sources  received  is  A.  Let  the  origin  of  a  coordinate  system  be  at  the 
mid-point  of  the  array,  with  the  sensors  individually  spaced  by  d  regularly  along  the  x-axis 
and  indexed  i  =  1, . . . ,  M  from  left  to  right.  We  assume  that  d<  f-  Boresight  is  along 
the  y-axis. 

A  source  is  placed  in  the  near-field  (i.e.  a  circular  wavefront  impinges  on  the  array) 
at  location  rs  meters  from  the  origin  and  6S  degrees  from  boresight.  For  convenience 
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(although  somewhat  unconventionally)  we  have  defined  that  angles  are  to  be  measured 
clockwise  from  array  boresight  (the  y-axis).  For  M  odd  there  is  a  sensor  at  the  origin, 
whereas  for  M  even  the  origin  is  midpoint  between  two  sensors.  The  array  geometry  and 
the  notations  are  shown  in  Fig.  10  for  the  case  of  M  even. 


y 


source 


Fig.  10.  Sensor  array  geometry  and  notations  for  linear  array  and  a  near-field  source 


The  distance  from  the  ith  sensor  to  the  source  is  given  by 


ri  +(*-!)• d 


M-  1 


•  d\  —  2rs  •  d-  (t  -  1) 


M-  1 


sin(0s)  (1) 


and  the  corresponding  complex  response  at  the  ith  sensor  is 


(  1  /  .27 r 

ai{rs,es )  = - exp  (  — j  —  *  rSi, 

rs.i  V  A 


assuming  a  normalized  and  equal  gain  for  each  sensor.  The  vector  a(rs,  6S)  =  [ai, . . . ,  aM]T 
is  the  response  of  the  complete  array  to  the  source  at  (rs,  0S). 

Likewise,  the  angle  from  the  ith  sensor  to  the  source  is 


6S  i  =  cos' 


1  TC(» -  !)  —  Mrf^  +  rh-rr\  tt 


2d[{i  - 1)  - 


Mnilr  • 

2  Jfs.1 


It  is  a  characteristic  of  near-field  sources  that  they  are  viewed  at  different  angles  by  the 
different  sensors  in  the  array. 

Given  sensor-angle  measurements  from  any  two  or  more  sensors,  it  is  possible  to  deter¬ 
mine  the  range  and  bearing  (rs,0s)  (with  respect  to  the  origin)  of  a  source  in  the  near-field. 
This  is,  however,  subject  to  identifiability  requirements  that  each  sensor  has  a  different 
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angle  to  the  source.  Given  0Sji  and  9sj  with  i  /  j,  one  determines  sensor  to  source  ranges 
rS)j  and  rsj  respectively  using 

(4) 

and 

This  requires  that  6s>i  —  0sj  /  0.  The  range  and  bearing  with  respect  to  the  origin  can  be 
determined  relative  to  any  of  the  individual  sensors  using  the  individual  sensor  range  and 
bearing.  For  example,  for  the  jth  sensor,  we  use  rsj  and  9sj  according  to 


_  2  i 

rs  ~~  rs,j  d" 


r  M-l  1“  rjw-l  1 

l^Y1  -  - 2  •  •  [[^y1  ■ [j  - 1]]  •  d\ sin[0 


M-l 

2 


0S  =  —  cos-1[0sJ].  (7) 

rs 

A. 2  Model 

Our  proposed  source  characterization  technique  requires  one  cooperative  source  with 
complex  envelope  sf  in  the  far-field  of  the  array  at  known  angle  9s.  Steering  vectors  for 
the  far-field  a(0)  and  near-field  a(0,  r)  take  on  the  standard  form  with  0  being  the  angle, 
r  denotes  range  [29]. 

Assume  that  the  conditions  on  the  test  source  and  sensor  array  are  such  that  the  fol¬ 
lowing  signal  model  is  appropriate 

zk  =  Ask  +  qk  +  nk.  (8) 

In  this  model,  zk  is  the  kth  snapshot  of  sensor  data  outputs  (dimension  M).  The  qk 
represents  additive  spatial  and  temporal  colored  noise  produced  in  the  environment  and 
nk  represents  additive  white  noise  modeling  the  internal  noise  of  the  array  of  sensors 
receiving  system. 

The  matrix  A  can  take  on  one  of  two  forms,  depending  on  whether  the  local  scatterer 
is  best  modeled  as  a  collection  of  P  discrete  scatterers,  or  as  a  single  spatially  distributed 
scatterer.  For  the  case  of  P  discrete  scatterers 


A  =  ja(0s), a(#i, Ti), . . . , a(0P, rP)] . 
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In  the  above  equation,  the  near-field  scatterer  i  =  1, . . . ,  P  is  characterized  by  the  angle 
6i  and  range  r*.  For  a  distributed  scatterer  with  scatter  amplitude  /[a(0,  r)]  contained  in 
the  near-field  azimuth  and  range  set  Q, 

A=[a (6s),  [  f  [a(0,  r)]d0drl .  (10) 

Likewise,  the  signal  vector  Sk  can  be  constructed  in  two  ways,  depending  on  whether 
the  near-field  scatterers  are  best  modeled  as  discrete  or  continuous.  For  the  discrete  case, 

Sk=[ssk,sL...,skp]T  (11) 

where  the  test  source  complex  amplitude  is  given  by  sf  and  the  complex  amplitude  of  the 
jtk  0£  p  discrete  scatterers  is  denoted  as  s\.  In  the  continuous  scatterer  case, 

Sk  =  [sk,Sk]T-  (12) 

The  4  and  slk  and  sk  may  be  uncorrelated,  correlated,  or  coherent  for  each  case  respec¬ 
tively. 

The  spatial  covariance  matrix  E[zzH]  is 

R  =  ASAh  +  Q  +  <j2I.  (13) 


Let  the  elements  of  S  be  p.tJ .  We  ensure  that  the  cooperative  test  source  has  sufficient  signal 
to  noise  ratio  (generally  greater  than  50dB)  to  perform  our  measurement  by  requiring  that 


Psnr  — 


(a2  +  tr[Q]) 


»  1. 


It  is  also  expected  that  the  direct  far-field  source  power  will  be  substantially  greater  than 
the  total  near-field  power  (by  20-40dB) 


A. 3  Background 


Psn(  = 


Pn 

P 

2  Pk 


»  1. 


(15) 


Breed  and  Posch  [27]  introduced  the  spatial  Wigner  distribution  as  a  tool  for  deter¬ 
mining  the  range  and  angle  of  a  near-field  source.  They  exploited  the  property  that  the 
phase  front  of  a  wave  emanating  from  a  source  in  the  near-field  of  an  array  has  an  approx¬ 
imately  quadratic  phase  law,  or  equivalently  an  approximately  linear  spatial  frequency 
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law.  They  then  determined  source  location  by  determining  the  parameters  governing  the 
linear  frequency  law  as  represented  by  the  Wigner  distribution  applied  to  the  spatial  sig¬ 
nal.  The  true  propagating  wave  phase  front  is  in  fact  spherical  and  is  only  approximately 
quadratic  for  near-field  sources  some  distance  from  the  array.  The  method  proposed  in  [27] 
breaks  down  for  sources  close  to  the  array.  Swindlehurst  and  Kailath  [30]  examined  the 
applicability  of  the  quadratic-spherical  approximation  and  apply  a  parametric  high  resolu¬ 
tion  technique  to  determine  the  linear  frequency  law  parameters  (and  hence  the  near-field 
source  position).  However,  it  can  be  seen  from  equations  (4)  through  (7)  that  it  is  possible 
to  determine  the  source  position  without  invoking  the  quadratic  phase  approximation  to 
the  spherical  phase  front. 

There  is  a  substantial  body  of  literature  concerned  with  processing  spatial  signals  re¬ 
ceived  by  an  array  of  sensors  from  sources  in  the  near-field  of  the  array.  It  is  mostly 
concerned  with  techniques  for  estimating  the  angle  and  range  of  the  source.  For  example, 
both  subspace  and  maximum  likelihood  algorithms  are  derived  in  [31].  Subsequently  we 
will  present  an  example  showing  near-field  characterization  using  both  the  sensor-angle 
distribution  (discussed  next)  and  the  near-field  MUSIC,  as  developed  in  [31]. 

Several  authors  have  proposed  methods  for  determining  the  angle  of  distributed  sources 
located  in  the  far-field  of  an  array  [32].  These  techniques  address  the  effect  of  scatter  local 
to  a  transmitter  in  the  far-field  and  not  for  scatter  that  is  sufficiently  local  to  the  receiving 
system  to  be  in  the  near-field  of  the  array. 

B.  Sensor-Angle  Distributions 

Our  method  extends  the  spatial  Wigner  distribution  introduced  by  Breed  and  Posch.  To 
avoid  confusion  it  has  been  necessary  to  change  the  name  to  reflect  the  generalization  to 
all  members  of  Cohen’s  class  of  quadratic  distributions  [1].  While  the  title  “spatial  Wigner 
distribution”  is  informative,  retaining  the  name  “spatial  time-frequency  distribution”  for 
the  remaining  members  of  Cohen’s  class  applied  to  spatial  signals  does  not  correctly  de¬ 
scribe  the  distribution  we  are  interested  in,  and  will  be  confused  with  STFD  discussed 
in  earlier  sections  of  this  chapter.  Therefore,  in  this  work  we  have  renamed  the  class  of 
quadratic  distributions  applied  to  spatial  signals  to  be  sensor-angle  distributions  (SAD). 
The  corresponding  spectra  are  called  sensor-angle  spectra  (SAS). 
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The  Cohen’s  class  SAD  for  the  kth  snapshot  is  a  distribution  of  the  angle  of  sources 
impinging  on  the  array  at  each  sensor. 

OO  OO 

Tk(i,  9;  zk)  =  £  <^(v,l)zk(i  +  v  +  l)Zk(i  +  v-l)e_j4,rM  (16) 

V—  oo  l=-oo 

where  %  and  6  are  the  sensor  index  and  angle  respectively.  The  kernel  <j)(v,  l )  characterizes 
the  distribution  and  is  a  function  of  sensor  position  and  sensor  lag.  All  the  standard  kernel 
designs  applied  in  the  time-frequency  literature  may  be  used  with  the  SAD. 

The  sensor-angle  spectrum  (SAS)  is  the  power  (not  energy  or  energy  density)  distribu¬ 
tion  of  the  sources  impinging  on  the  array.  The  SAS  is  given  by 

Ts(i,M  =  E[Tk(i,0;i!k)]  (17) 

where  an  estimate  for  temporally  stationary  sources  is  given  by 

fs(i,0;z)  =  l^1Tk(i,0;zk)  (18) 

iN  k=0 

for  N  snapshots. 

C.  Characterizing  Local  Scatter 

The  objective  is  to  use  data  received  by  the  array  from  a  test  source  in  the  far-held  that 
illuminates  the  local  near-held  scatterer  distribution  and  to  visualize  and  characterize  this 
scatterer  distribution  using  the  SAS.  We  expect  the  test  signal  to  be  substantially  more 
powerful  than  the  local  scatter  we  wish  to  characterize  (see  (15)).  Subspace  projection  is 
applied  to  the  array  snapshots  to  remove  the  dominant  far-held  component  and  allow  a 
clear  depiction  of  the  near-held  source  in  the  sensor-angle  (s-a)  domain. 

In  (16)  and  (17),  the  data  snaphot  zk  is  replaced  by  PflSzk  where  PeS  is  the  orthogonal 
projection  operator  formed  from  the  far-held  test  source  steering  vector  a (0s)  as 

P*S  =  I  -  a(0s)[aH(0s)a(0s)]-1aH(0s).  (19) 

Therefore,  we  compute  the  modihed  SAS 

i^MiP*8,*).  (20) 

In  some  applications,  a  single  test  angle  will  provide  sufficient  characterization  using 
(20)  while  in  other  applications,  two  or  several  test  angles  will  be  required,  in  which  case 
0s  is  scanned  over  the  required  domain  of  angles  for  the  test  source. 


84 


D.  Simulations  and  Examples 

The  following  example  is  used  to  demonstrate  the  proposed  approach  for  near  scattering 
characterizations.  Consider  a  32  sensor  linear  equi-spaced  array  operating  at  a  carrier 
frequency  of  6.41MHz  and  with  15m  sensor  spacing.  The  local  scatterer  distribution 
comprises  a  point  scatterer  in  the  near-field  at  a  range  of  400m  and  bearing  of  30  degrees 
in  front  of  the  array  (the  array  has  total  length  465m).  Assume  the  test  source  is  temporally 
stationary  and  located  at  20  degrees  angle  with  respect  to  boresight.  The  test  source  is 
coherent  with  and  20dB  stronger  than  the  scattered  source.  In  this  example  we  have  used 
the  alias-free  Wigner  distribution  [33].  Of  course  others  members  of  Cohen’s  class  may 
also  be  used. 

Figure  11  shows  the  SAD  for  the  received  data.  The  SAD  is  dominated  by  the  substan¬ 
tially  more  powerful  far-field  test  source  and  there  is  no  clear  indication  of  any  additional 
scattering.  The  far-field  source  has  the  same  angle  for  every  sensor,  and  therefore,  depicts 

a  horizontal  signature  in  the  s-a  domain.  In  Fig.  12  we  have  applied  the  orthogonal  projec- 

s 

tion  operator  and  computed  the  SAD  for  P6  zk.  The  SAD  now  clearly  shows  the  presence 
of  near-field  local  scatter.  The  location  of  the  near-field  source  may  be  determined  using 
equations  (4)-(7). 

Sansor-angto  distribution  (alias-free  Wigner  dist) 


I 

! 

1 


Fig.  11.  SAD  for  the  received  data  zk.  The  far-field  test  source  dominates  the  SAD  characterization. 

The  beampatterns  for  the  cases  of  zk  and  P^Szk  are  shown  in  Fig.  13.  The  presence 
of  near-field  scatterers  cannot  be  confirmed  as  compared  with  alternative  explanations  for 
the  distorted  beampatterns,  such  as  poor  array  calibration. 
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Fig.  12.  SAD  for  the  received  data  P*  zk.  With  the  direct  propagation  path  from  the  far-field  test  signal 
removed  by  the  orthogonal  projection  operator  the  local  scatterer  spectrum  is  revealed. 

Baampattam 


g 

Fig.  13.  Beampattern  for  Zk  ( — )  and  P°  Zk  ( — ).  Without  the  sensor-angle  characterization  it  is  not 
possible  to  identify  perturbations  from  the  ideal  test  source  beampattern  as  being  due  to  near-field 
scatter. 


The  projection  approach  has  also  been  applied  to  real  data  collected  from  a  16  sensor 
HF  receiving  array.  An  array  calibration  source  was  transmitted  from  the  far-field  of  the 
array  at  boresight.  The  sensor  angle  distribution  is  shown  in  Fig.  14  and  is  dominated 
by  the  calibration  source.  Following  calibration  of  the  array  using  the  calibration  source, 
the  received  and  calibrated  boresight  source  is  removed  using  orthogonal  projection.  The 
sensor  angle  distribution  of  the  residual  is  shown  in  Fig.  15.  No  discrete  near-field  scatterers 
are  apparent,  however  there  is  a  concentration  in  the  SAD  in  the  upper  left  region  of  the 
distribution.  This  indicates  that  there  is  some  asymmetric  local  scattering  near  the  array. 

A  second  example  is  used  to  contrast  the  SAD  with  existing  techniques  for  near-field 
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Fig.  14.  SAD  for  the  real  received  data.  The  far-field  calibration  source  at  boresight  dominates  the  SAD 
characterization. 


Sensor-angle  distribution  (alas— free  Wigner  dist) 
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Fig.  15.  SAD  for  the  real  received  data  with  the  direct  propagation  path  from  the  far-field  calibration 
source  removed  by  orthogonal  projection  operator.  The  local  scatterer  sensor-angle  distribution  is 
revealed. 

sources  characterization.  In  this  case,  we  have  chosen  to  compare  with  an  implementation 
of  near-field  MUSIC  as  described  by  [31].  We  consider  an  ideal  computer  generated  case 
and  an  equivalent  case  where  the  data  has  been  collected  using  a  real  HF  radar  array. 

Consider  an  ideal  case  with  a  single  source  placed  at  range  80m  and  bearing  10.5  degrees. 
Assume  a  16  sensor  array,  that  source  signal  to  noise  power  ratio  is  high  and  102  array  data 
snapshots  are  available  (generated  using  computer).  Figs.  16  and  17  respectively,  show 
the  near-field  MUSIC  diagram  and  the  SAS,  both  computed  using  all  102  data  snapshots. 
The  source  is  well  localized  using  MUSIC,  and  there  is  a  characteristic  structure  in  the 
SAS  showing  sensor-angle  for  each  sensor  in  the  array. 
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Nsar-fiaM  MUSIC  orthogonality  measure  (dB  wrt  max) 


angle -ot-a»rival  (deg) 

Fig.  16.  Near-field  MUSIC  diagram  for  an  ideal  source  at  range  80m  and  bearing  10.5  degrees. 


Sensor-angle  spectrum  (alias -free  Wtgner  tfct)  1 02  snapshots 
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Fig.  17.  sensor-angle  spectrum  for  an  ideal  source  at  range  80m  and  bearing  10.5  degrees. 

We  have  repeated  the  analysis,  but  this  time  used  102  data  snapshots  collected  from  the 
HF  radar  receiving  array.  The  near-field  source  was  approximately  80m  from  the  array 
mid-point  at  an  angle  of  approximately  10  degrees.  In  this  case  the  source  was  behind 
the  array.  The  exact  location  is  not  known  precisely.  Figs.  18  and  19  show  the  near-field 
MUSIC  diagram  and  the  SAS  respectively.  Imprecise  array  calibration  has  smeared  the 
localization  in  the  MUSIC  diagram  while  the  structure  in  the  SAS  is  preserved. 


VIII.  Conclusion 

We  have  presented  two  different  new  perspectives  of  time-frequency  distributions.  One 
perspective  is  driven  by  direction  finding  and  blind  source  separation  problems,  whereas 
the  other  stems  from  the  need  to  characterize  near-field  sources  or  reflectors.  The  fun- 


N«ar -field  MUSIC  orthogonality  measure  (dB  wrtmax) 


-80  -60  -40  -20  0  20  40  60  80 

angle-of-airivat  (deg) 


Fig.  18.  Near-field  MUSIC  diagram  for  a  real  source  placed  at  approximately  range  80m  and  bearing 
10.5  degrees. 


Sensor-angle  spectrum  (alias-free  Wigner  dbt)  102  snapshots 
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Fig.  19.  sensor-angle  spectrum  for  a  real  source  placed  at  approximately  range  80m  and  bearing  10.5 
degrees. 

damental  offering  of  quadratic  distributions  in  both  cases  is  the  ability  to  discriminate 
between  the  sources  based  on  the  joint-variable  signatures  of  their  respective  waveforms. 
This  allows  the  enhancement  of  signal-to-noise  ratio  (SNR)  as  well  as  the  consideration 
of  only  the  sources  of  interest,  and  subsequently  improve  the  estimation  of  the  source 
positions  and  waveforms. 

The  first  six  sections  presented  the  general  framework  of  spatial  time-frequency  distri¬ 
butions  (STFDs).  The  advantages  of  a  STFD  matrix  over  the  covariance  matrix-based 
approach  to  array  processing  are  the  SNR  enhancement  and  the  robustness  of  the  eigen- 
structure  to  noise.  A  variety  of  methods  have  been  introduced  for  both  blind  source 
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separation  and  direction  finding  applications  using  STFDs.  The  first  class  of  source  sepa¬ 
ration  methods  are  based  on  pre- whitening  and  the  recovery  of  a  unitary  matrix.  Unlike 
similar  methods  based  on  second-order  statistics,  which  cannot  separate  signals  with  the 
same  spectra,  the  STFD-based  method  can  separate  nonstationary  signals  with  identical 
spectra  when  they  have  different  time-frequency  signatures.  The  SNR  enhancement  and 
signal  localization  properties  in  the  time-frequency  domain  can  substantially  improve  the 
source  separation  performance.  The  second  class  of  source  separation  methods  are  based 
on  array  averaging  of  the  time-frequency  distribution,  signal  synthesis,  and  waveform  re¬ 
covery  using  the  minimum  mean  square  error  criterion.  For  direction  finding,  both  the 
MUSIC  and  the  maximum  likelihood  methods  have  been  extended  and  modified  to  incor¬ 
porate  the  STFDs.  The  performance  improvement,  evident  by  both  the  analytical  and 
simulation  results,  is  most  significant  when  the  input  SNR  of  source  arrivals  is  low,  and/or 
when  the  sources  are  closely  spaced. 

In  the  second  part  of  the  chapter,  we  have  used  the  time-frequency  distribution  of  the 
spatial  signal  received  by  an  array  to  characterize  sources  based  on  their  angle  at  each 
sensor.  This  sensor-angle  distribution  is  a  tool  for  characterizing  near-field  scatter  local 
to  the  receiving  array.  The  method  uses  a  test  source  in  the  far-field  to  illuminate  the 
local  scatterer  distribution.  An  orthogonal  projection  operator  derived  from  the  steering 
vector  for  the  far-field  test  source  is  used  to  exclude  the  direct  propagation  path  from 
the  test  source  in  the  characterization.  As  part  of  the  characterization  we  exploit  the 
spatial  Wigner  distribution  although  we  have  renamed  it  the  sensor-angle  distribution  to 
avoid  confusion  with  a  similarly  named  but  differently  defined  STFD  discussed  in  the  first 
part  of  the  chapter.  We  have  shown  the  application  of  the  method  using  simulation  and 
for  real  data  collected  using  an  HF  radar  receiving  array.  Additional  simulations  and 
real  data  results  contrast  the  SAD  characterization  with  that  of  a  conventional  near-field 
localization  technique  (in  this  case  near-field  MUSIC). 
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Abstract 

This  chapter  discusses  interference  excision  techniques  in  spread  spectrum  communi¬ 
cations  systems.  Both  stationary  and  nonstationary  interferers  are  considered.  Sinusoidal, 
auto-regressive,  digital  communication,  and  polynomial  phase  interference  signals  are  ef¬ 
fectively  suppressed  by  different  excision  methods;  each  is  suitable  for  one  or  more  types 
of  interferers.  Time-,  frequency-,  and  time-frequency  domain  methods  for  analysis  and 
estimation  of  the  interference  parameters  are  summarized.  Domains  other  than  time  and 
frequency,  such  as  the  Gabor-domain,  the  Wavelet-domain  ,  and  quadratic  time-frequency 
signal  representations,  are  appropriate  for  non-traditional  smart  jamming  in  which  the 
interference  parameters  are  highly  dependent  on  time.  Excision  methods  can  be  linear, 
bilinear,  or  nonlinear  with  performance  dependent  on  the  interference  power  relative  to 
the  desired  signal  and  noise.  In  some  cases,  the  slight  performance  enhancement  offered 
by  bilinear  and  nonlinear  methods  may  not  properly  justify  the  increase  in  the  algorithm 
complexity.  The  BER  and  receiver  SINR  expressions  and  curves  are  presented  in  this 
chapter  for  some  of  the  key  interference  excision  techniques. 
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I.  Introduction 


Suppression  of  correlated  interference  is  an  important  aspect  of  modern  broadband  com¬ 
munication  platforms.  For  wireless  communications,  in  addition  to  the  presence  of  benign 
interferers,  relatively  narrowband  cellular  systems,  employing  time  division  multiple  ac¬ 
cess  (TDMA)  or  advanced  mobile  phone  system  (AMPS)  may  coexist  within  the  same 
frequency  band  of  the  broadband  code  division  multiple  access  (CDMA)  systems.  Hostile 
jamming  is  certainly  a  significant  issue  in  military  communication  systems.  Global  posi¬ 
tioning  system  (GPS)  receivers  potentially  experience  a  mixture  of  both  narrowband  and 
wideband  interference,  both  intentionally  and  unintentionally  [1]. 

One  of  the  fundamental  applications  of  spread  spectrum  (SS)  communication  systems  is 
its  inherent  capability  of  interference  suppression.  SS  systems  are  implicitly  able  to  provide 
a  certain  degree  of  protection  against  intentional  or  unintentional  interferers.  However, 
in  some  cases,  the  interference  might  be  much  stronger  than  the  SS  signal,  and  the  limi¬ 
tations  on  the  spectrum  bandwidth  render  the  processing  gain  insufficient  to  decode  the 
useful  signal  reliably.  For  this  reason,  signal  processing  techniques  are  frequently  used 
in  conjunction  with  the  SS  receiver  to  augment  the  processing  gain,  permitting  greater 
interference  protection  without  an  increase  in  the  bandwidth.  Although  much  of  the  work 
in  this  area  has  been  motivated  by  the  applications  of  SS  as  an  anti-jamming  method  in 
military  communications,  it  is  equally  applicable  in  commercial  communication  systems 
where  SS  systems  and  narrowband  communication  systems  may  share  the  same  frequency 
bands. 

This  article  covers  both  the  direct-sequence  spread  spectrum  (DS/SS)  and  frequency 
hopping  (FH)  communication  systems,  but  the  main  focus  is  on  the  DS/SS  communica¬ 
tion  systems.  For  DS/SS  communication  systems,  two  types  of  interference  signals  are 
considered,  namely,  narrowband  interference  (NBI)  and  nonstationary  interference,  such 
as  instantanously  narrowband  interference  (INBI). 

The  early  work  on  narrowband  interference  rejection  techniques  in  spread  spectrum 
communications  has  been  reviewed  comprehensively  by  Milstein  in  [2].  Milstein  discusses 
in  depth  two  classes  of  rejection  schemes:  1)  those  based  on  least-mean  square  (LMS) 
estimation  techniques,  and  2)  those  based  on  transform  domain  processing  structures.  The 
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improvement  achieved  by  these  techniques  is  subject  to  the  constraint  that  the  interference 
be  relatively  narrowband  with  respect  to  the  DS/SS  signal  waveform.  Poor  and  Rusch  [3] 
have  given  an  overview  of  NBI  suppression  in  SS  with  the  focus  on  CDMA  communications. 
They  categorize  CDMA  interference  suppression  by  linear  techniques,  nonlinear  estimation 
techniques,  and  multiuser  detection  techniques  (multiuser  detection  is  outside  the  scope 
of  this  article).  Laster  and  Reed  [4]  have  provided  a  comprehensive  review  of  interference 
rejection  techniques  in  digital  wireless  communications,  with  the  focus  on  advances  not 
covered  by  the  previous  review  articles. 

Interference  suppression  techniques  for  nonstationary  signals,  such  as  INBI,  have  been 
summarized  by  Amin  and  Akansu  [5].  The  ideas  behind  NBI  suppression  techniques  can 
be  extended  to  account  for  the  nonstationary  nature  of  the  interference.  For  time-domain 
processing,  time-varying  notch  filters  and  subspace  projection  techniques  can  be  used  to 
mitigate  interferes  characterized  by  their  instantaneous  frequencies  and  instantaneous 
bandwidths.  Interference  suppression  is  achieved  using  linear  and  bilinear  transforms, 
where  the  time-frequency  domain  and  wavelet/Gabor  domain  are  typically  considered. 
Several  methods  are  available  to  synthesize  the  nonstationary  interference  waveform  from 
the  time-frequency  domain  and  subtract  it  from  the  received  signal. 

Interference  rejection  for  FH  is  not  as  well  developed  as  interference  rejection  for  DS  or 
for  CDMA.  In  FH  systems,  the  fast  FH  (FFH)  is  of  most  interest,  and  the  modulation  most 
commonly  used  in  FH  is  frequency-shift  keying  (FSK).  Two  types  of  interference  waveforms 
can  be  categorized,  namely,  partial-band  interference  (PBI)  and  multitone  interference 
(MTI).  Typically,  interference  suppression  techniques  for  FH  communication  systems  often 
employ  a  whitening  or  clipping  stage  to  reject  interference,  and  then  combined  by  diversity 
techniques. 

II.  Signal  Model 

The  received  waveform  consists  of  a  binary  phase-shift- keying  (BPSK)  DS/SS  signal 
s(t),  an  interfering  signal  u(t),  and  thermal  noise  b(t).  Without  loss  of  generality,  we 
consider  the  single-interferer  case,  and  additive  Gaussian  white  noise  (AGWN)  that  is 
uncorrelated  with  both  the  DS/SS  and  the  interference  signals.  The  input  to  the  receiver, 
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x(t),  is  given  by 

x(t)  —  s(t)  +  u(t )  +  b(t).  (1) 

The  DS/SS  signal  can  be  expressed  as 

OO 

s(t)  =  ■£  d(l) p(t  -  ITC),  (2) 

/  —  —  OO 

where  Tc  is  the  chip  duration,  p(t)  is  the  shaping  waveform, 

d(l)  =  s(ri)c(n,  l )  (3) 

is  the  chip  sequence,  s(n)  G  [—1,  +1]  is  the  nth  symbol,  and  c(n,  l)  G  [—1,  +1]  is  a  pseudo¬ 
noise  (PN)  sequence  used  as  the  spreading  code  for  the  nth  symbol.  The  PN  sequence  can 
be  either  periodic  and  aperiodic.  Different  types  of  interference  signals  are  considered. 

For  discrete-time  filter  implementations,  signals  are  sampled  at  the  rate  1/T.  Typically, 
the  sampling  interval  T  is  equal  to  the  chip  duration  Tc.  The  input  to  the  receiver,  after 
sampling,  becomes 

x[n]  =  x(nT).  (4) 

The  samples  of  the  DS/SS  signal,  interference,  and  noise,  can  be  defined  accordingly  as 
s[n],  «[n],  and  b[n],  respectively. 

III.  Narrowband  Interference  Suppression 

Interference  suppression  techniques  for  DS/SS  systems  are  numerous.  In  particular, 
much  literature  exists  on  the  adaptive  notch  filtering  as  it  relates  to  suppress  NBI  on  a 
wideband  DS/SS  signal.  Synthesis/subtraction  is  another  well-established  technique  for 
sinusoidal  interference  suppression.  Other  techniques  include  nonlinear  adaptive  filtering 
and  multiuser  detection  techniques. 

A.  Adaptive  Notch  Filtering 

The  basic  idea  in  employing  an  adaptive  notch  filter  is  to  flatten  the  filter  input  spec¬ 
trum.  An  SS  signal  tends  to  have  a  uniform  wide  spectrum  and  is  affected  little  by  the 
filtering  process,  whereas  the  NBI  is  characterized  by  spectral  spikes  and  frequency  regions 
of  high  concentrated  power.  The  adaptive  notch  filter  places  notches  at  the  frequency  lo¬ 
cation  of  the  NBI  to  bring  the  interference  level  to  the  level  of  the  SS  signal.  At  least  two 
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main  approaches  exist  for  constructing  an  adaptive  notch  filter:  (1)  estimation/subtraction 
type  filters  and  (2)  transform-domain  processing  structures. 

Estimation/Subtraction  Type  Filters 

If  the  interference  is  relatively  narrowband  compared  with  the  bandwidth  of  the  spread 
spectrum  waveform,  then  the  technique  of  interference  cancellation  by  the  use  of  notch 
filters  often  results  in  a  large  improvement  in  system  performance.  This  technique,  de¬ 
scribed  in  many  references  including  [6],  [7],  [8],  [9],  uses  a  tapped  delay  line  to  implement 
the  prediction-error  filter  (Wiener  filter  [10]).  Since  both  the  DS  signal  and  the  thermal 
noise  are  wideband  processes,  their  future  values  cannot  be  readily  predicted  from  their 
past  values.  On  the  other  hand,  the  interference,  being  a  narrow-band  process,  can  indeed 
have  its  current  and  future  values  predicted  from  past  values.  Hence,  the  current  value, 
once  predicted,  can  be  subtracted  from  the  incoming  signal,  leaving  an  interference-free 
waveform  comprised  primarily  of  the  DS  signal  and  the  thermal  noise.  A  general  diagram 
of  this  technique  is  depicted  in  Fig.  1.  Both  one-sided  and  two-sided  transversal  filters 
can  be  used  for  this  purpose.  When  two-sided  filters  are  used,  the  estimation  of  current 
interference  value  is  based  on  both  past  and  future  values  of  the  interference.  Consider  a 
single-sided  filter  as  shown  in  Fig.  2.  Define  an  iV-dimensional  vector  x[n],  denoted  as 

x[n]  =  (x[n  —  1],  •  •  • ,  x[n  —  N])T  ,  (5) 

where  the  superscript  T  denotes  transpose  of  a  vector  or  a  matrix.  The  DS/SS  signal, 
interference,  and  noise  vectors  can  be  defined  similarly  as  s[n],  u[n],  b[n],  respectively.  We 
also  define  the  corresponding  weight  vector  w  as 

w  =  [uq,-  ••,w/v]T.  (6) 

Hence,  the  output  sample  of  the  filter  is 

y[n]  =  x[n]  —  wTx[n].  (7) 

The  mean  square  value  E[y2[n]],  representing  the  output  average  power,  is  given  by 

E  (y2[n])  =  E  (x2[n])  -  2w TE  (x[n]x[n])  +  w TE  (x[n]xT[n])  w 

A  E  (x2[n])  -  2wTp  +  wtRw,  (8) 
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where  p  =  E  (x[n]x[n])  is  the  correlation  vector  between  x[n]  and  x[n],  and 

R  —  E  (x[n]xr[n])  (9) 

is  the  covariance  matrix  of  x[n].  It  is  noted  that,  when  the  PN  sequence  is  sufficiently 
long,  the  PN  signal  samples  at  different  taps  are  approximately  uncorrelated.  On  the 
other  hand,  samples  of  the  narrowband  interference  at  different  taps  has  high  correlations. 
Since  the  DS/SS  signal,  interference,  and  thermal  noise  are  mutually  uncorrelated,  then, 

p  =  E  (a:[n]x[n]) 

=  E  |(s[n]  +  u[n]  +  b[n ])  (s[n]  +  u[n]  +  b[n])} 

=  E  (u[n]u[n]) .  (10) 

Minimizing  the  output  power  E[y2[n ]]  yields  the  following  well-known  Wiener-Hopf  so¬ 
lution  for  the  optimum  weight  vector  wopt, 

wopt  =  R_1p.  (11) 


The  cost  of  notch  filtering  is  the  introduction  of  some  distortion  into  the  SS  signal.  Such 
distortion  results  in  power  loss  of  the  desired  DS/SS  signal  as  well  as  the  introduction  of 
self-noise.  Both  effects  become  negligible  when  the  processing  gain  is  sufficiently  large. 

Note  that  when  precise  statistical  knowledge  of  the  interference  cannot  be  assumed, 
adaptive  filtering  can  be  used  to  update  the  tap  weights.  There  are  a  variety  of  adaptive 
algorithms  and  receiver  structures  [7],  [11],  [12],  [10].  The  optimum  Wiener-Hopf  filter  can 
be  implemented  by  using  direct  matrix  inversion  (DMI)  or  recursive  adaptation  methods. 
For  the  DMI  method,  the  covariance  matrix  R  and  p  are  estimated  at  time  n  by  using 
most  recent  Nt  data  samples,  i.e., 

1  Nt- 1 

R-W  =  T7  ]£  Xtn  “  ^[n  ~  *1  (12) 

1=0 


1  N‘~l 

PM  =  TF  S  4n  ~  -  1 ]  (13) 

1=0 

The  least  mean  square  (LMS)  algorithm  is  a  simple  and  stable  method  to  implement 
an  iterative  solution  to  the  Wiener-Hopf  equation  without  making  use  of  any  a  priori 
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statistical  information  about  the  received  signal.  Using  the  instantaneous  estimates  of  the 
covariance  matrix  and  cross-correlation  vector  of  Eqns.  (9)  and  (10),  the  LMS  algorithm 
can  be  expressed  as 

w[n  +  1]  =  w[n]  +  /xy[n]x[n],  (14) 

where  w[n]  is  the  filter  weight  vector  of  elements  W[ [n],  l  —  1,2,*--,  AT,  x[n]  is  the  vector 
that  includes  the  data  within  the  filter,  and  y[n]  is  the  output,  all  at  the  nth  adaptation, 
and  n  is  a  parameter  that  determines  the  rate  of  convergence  of  the  algorithm.  It  is  noted 
that  in  most  applications  of  the  LMS  algorithm,  an  external  reference  waveform  is  needed 
in  order  to  correctly  adjust  the  tap  weights.  However,  in  this  particular  application,  the 
signal  on  the  reference  tap  x[n]  serves  as  the  external  reference. 

The  drawback  of  LMS  algorithm  is  its  slow  convergence.  To  improve  the  convergence 
performance,  techniques  including  self-orthogonalizing  LMS,  recursive  least  squares  (RLS), 
and  lattice  structure  can  be  used  [13],  [10]. 

SINR  and  BER  Analysis 

The  output  signal-to-interference-plus-noise  ratio  (SINR)  and  bit  error  rate  (BER)  are 
two  important  measures  for  communication  quality  and  the  performance  enhancement 
using  the  signal  processing  techniques. 

To  derive  the  output  SINR,  we  rewrite  the  filter  output  as  1 

N  N 

y[n  ]  =  'jT  wix[n  —  /]  =  ^  wi{c[n  —  /]  +  u[n  —  /]  +  b[n  —  I]).  (15) 

1=0  1=0 

where  w0  =  1.  The  signal  {y[n]}  is  then  fed  to  the  PN  correlator.  Denote  L  as  the  number 
of  chips  per  information  bit.  Then  the  output  of  the  PN  correlator,  which  is  the  decision 
variable  for  recovering  the  binary  information,  is  expressed  as 

LLN 

r  —^2  y[n]c[n }  -  ^  c[n ]  ^  wi(c[n  —  I]  +  u[n  -  /]  +  b[n  -  /]) 

n= 1  7i—l  /=0 

=  £  c2M  +  5Z  c[n]  X]  wicin  +  cM  53  wiiu\n  - 1]  +  bin  -  *]) 

n=  1  7i—l  l— 1  Ti—1  /— 0 

L  N  L  N  L  N 

=  L  +  ^2  c[n]  ^2  wtc[n  -  /]  +  ^  X]  c[n]wiu[n  ~  l]  +  CW wMn  ~  A- 

n—l  1=1  7i=l  1=0  7i=l  1=0 

lrro  keep  the  notation  simple,  we  have  used  in  (15)  the  same  symbol  wi  as  in  (6).  The  two  sets  of  weights, 
however,  differ  in  sign. 
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(16) 


The  first  term  on  the  right-hand  side  of  (16)  represents  the  desired  signal  component,  the 
second  term  amounts  to  the  self-noise  caused  by  the  dispersive  characteristic  of  the  filter, 
and  the  third  term  is  the  residual  narrowband  interference  escaping  the  excision  process 
and  appearing  at  the  output  of  the  PN  correlator.  The  last  term  in  (16)  is  the  additive 
noise  component. 

The  mean  value  of  r  is 

E(r)  =  L  (17) 


and  the  variance  is  [7] 

var(r)  Act  2  =  L  53  wf  +  L  £  wnwip[n  -/]  +  Lai  2Z  w 'f  >  (18) 

1=1  n=l  1=0  1=0 

where  a2  =  E(b2[n])  is  the  AGWN  variance  and 


p[n  — 1]  =  E  (u[k]u[k  +  n  —  /]) . 


The  three  terms  of  the  right-hand  side  of  (18)  represent  the  mean  square  values  caused  by 
the  self-noise,  residual  narrowband  interference,  and  noise,  respectively. 

The  output  SINR  is  defined  as  the  ratio  of  the  square  of  the  mean  to  the  variance.  Thus, 

SINR.  =  - ■  (19) 

1=1  n=l  1=0  1=0 

Note  that,  if  there  is  no  interference  suppression  filter,  wi  =  1  for  /  =  0  and  zero  otherwise. 
Therefore,  the  corresponding  output  SINR  is 

SINRn0  =  L  2.  (20) 

p[0]  +  ol 

If  we  assume  that  the  self-noise,  residual  interference,  and  noise  components  at  the 
output  of  correlator  is  Gaussian,  then  the  BER  can  be  evaluated  in  the  same  manner  as 
the  conventional  BPSK  corrupted  only  by  AWGN.  Under  such  assumption,  the  BER  is 
given  by 

Pb  =  P(r  <  0)  =  /°  -7L-e-(T-L)2/2°2dr  =  Q(\/Sl NRo)  ,  (21) 

7-oo  V27 xa  v  ' 
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where 


Qw=^rv’,/2ifo 


(22) 


is  the  Q-function  [13]. 

Transform- Domain  Processing  Structures 

An  alternative  to  time-domain  excision  as  described  in  the  preceding  section,  is  to  trans¬ 
form  the  received  signal  to  the  frequency  domain  and  perform  the  excision  in  that  domain. 
Clipping  and  gating  methods  can  then  be  applied  on  those  transform  bins  contaminated 
by  the  interference. 

Surface  acoustic  wave  (SAW)  device  technology  can  be  used  to  produce  the  continuous¬ 
time  Fourier  transform  of  the  received  waveform  [14],  [15].  The  discrete  Fourier  transform 
(DFT),  with  FFT  implementations,  is  commonly  applied  for  time-sampled  signals  [16]. 
Adaptive  subband  transforms  generalize  transform-domain  processing  [17],  [18],  and  can 
yield  uncorrelated  transform  coefficients. 

The  interference-suppressed  signal  based  on  a  block  transform  can  be  written  as 

xs[n]  =  BEAx[n],  (23) 

where  x[n]  is  the  received  input  vector,  A  and  B  are  the  forward  transform  matrix  and 
inverse  transform  matrix,  respectively,  and  E  is  a  diagonal  matrix  with  each  diagonal 
element  acting  as  a  weight  multiplied  to  the  input  signal  at  each  transform  bin.  The 
weights  can  be  controlled  by  different  schemes.  Two  commonly  used  methods  are  either 
to  set  the  weights  binary  (i.e.,  a  weight  is  either  one  or  zero)  or  to  adjust  the  weights 
adaptively.  In  applying  the  first  method,  powerful  NBI  is  detected  by  observing  the 
envelope  of  the  spectral  waveform.  Substantial  interference  suppression  can  be  achieved 
by  multiplying  the  input  signal  with  a  weight  that  is  set  to  zero  when  the  output  of  the 
envelope  detector  at  a  transform  bin  exceeds  a  predetermined  level.  Fig.  3  illustrates 
the  concept  of  transform-domain  notch  filtering.  Adaptive  algorithms,  such  as  LMS  and 
RLS,  can  be  used  to  determine  the  excision  weights  adaptively.  The  application  of  these 
algorithms,  however,  requires  a  reference  signal  that  is  correlated  with  the  DS/SS  signal. 

When  the  binary  weights  are  used,  the  transform-domain  processing  technique  may  suf¬ 
fer  from  the  interference  sidelobes.  With  block  transforms,  the  energy  in  the  narrowband 
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interference,  which  initially  occupies  a  small  region  of  the  frequency  spectrum,  is  dispersed 
in  a  relatively  large  spectral  region.  In  this  case,  excision  of  a  large  frequency  band  may 
become  necessary  to  effectively  remove  the  interference  power  from  most  transform  bins. 
The  frequency  dispersion  of  the  interference  can  be  nearly  eliminated  by  weighting  the 
received  signal  in  the  time  domain  with  a  nonrectangular  weighting  function  prior  to  eval¬ 
uating  the  transform.  In  doing  so,  the  levels  of  the  side  lobes  of  the  interference  frequency 
spectrum  are  attenuated  at  the  expense  of  broadening  the  main  lobe  [19],  [15].  In  this 
case,  the  conventional  matched  filter  is  no  longer  optimal.  Using  adapted  demodulation 
accordingly  can  improve  the  receiver  performance  [20]. 

It  is  important  to  point  out  that  for  transform-domain  processing,  symbol  detection  can 
be  performed  in  either  the  time  or  the  transform  domain.  In  the  later  case,  filtering  and 
correlation  operations  can  be  combined  in  one  step. 

The  BER  expression  for  transform-domain  interference  excision  can  be  easily  formulated 
using  the  Gaussian  tail  probability  or  the  Q-function.  The  residual  filtered  and  despreaded 
interference  is  treated  as  an  equivalent  AWGN  source.  Typically,  a  uniform  interference 
phase  distribution  is  assumed.  When  transform-domain  filtering  is  considered,  the  BER 
depends  on  both  the  excision  coefficients  and  the  error  misadjustment. 

B.  Synthesis/Subtraction  for  Sinusoidal  Interference 

In  this  section,  we  view  the  interference  signal  as  the  one  that  is  corrupted  by  the 
additive  noise  and  the  DS/SS  signal. 

Eigenanalysis  Approaches 

In  a  typical  situation,  the  power  level  of  the  DS/SS  signal  is  negligible  relative  to  the 
power  level  of  the  interference  and,  in  the  most  case,  relative  to  the  additive  noise.  For 
high  interference-to-noise  ratio  (INR),  the  correlation  matrix  of  the  received  signal  vector 
consists  of  a  limited  number  of  large  eigenvalues  contributed  mainly  by  the  narrowband 
interference,  and  a  large  number  of  small  and  almost  equal  eigenvalues  contributed  by 
the  DS/SS  signal  and  noise.  The  eigenanalysis  interference  canceller  is  designed  with  a 
weight  vector  orthogonal  to  the  eigenvectors  corresponding  to  the  large  eigenvalues  [21]. 
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The  eigendecomposition  of  the  correlation  matrix,  defined  in  (9),  results  in 

"S  0  1  rUff" 

E  =  USU«  =  [UrU„]  '  /  ,  (24) 

-  0  2jn  J  L  . 

where  the  columns  of  Ur  span  the  interference  subspace,  whereas  the  columns  of  Un  span 
the  signal  with  noise  subspace,  and  £r  and  £n  are  diagonal  matrices  whose  elements  are 
the  eigenvalues  of  R.  For  real  sinusoidal  interference,  the  number  of  dimensions  of  the 
interference  subspace  is  twice  the  number  of  interfering  tones. 

The  projection  of  the  signal  vector  on  the  noise  subspace  results  in  interference  sup¬ 
pressed  data  sequence, 

x[n]  =  UnU*x[n]  =  (i  -  UrUf )  x[n],  (25) 

where  I  is  the  identity  matrix. 

The  subspace  projection  approach  can  also  be  performed  using  the  singular  value  de¬ 
composition  (SVD)  for  the  sample  data  matrix  [22]. 

C.  Nonlinear  Estimation  Techniques 

The  commonly  applied  predictor/subtracter  technique  for  narrowband  interference  sup¬ 
pression  previously  discussed  is  optimum  in  the  minimum  mean  square  error  (MMSE) 
sense  when  trying  to  predict  a  Gaussian  autoregressive  process  in  the  presence  of  AWGN. 
If  the  prediction  is  done  in  a  non-Gaussian  environment,  as  in  the  case  of  SS  signals, 
linear  prediction  methods  will  no  longer  be  optimum.  In  [3],  depending  on  whether  the 
statistics  of  the  AR  process  is  known  or  unknown,  time-recursive  and  data-adaptive  non¬ 
linear  filters  with  soft  decision  feedback  are  used  to  estimate  the  SS  signal.  For  known 
interference  statistics,  the  interference  suppression  problem  is  cast  in  state  space  for  use 
with  Kalman-Bucy  and  approximate  conditional  mean  (ACM)  filters.  A  fixed  length  LMS 
transversal  filter,  on  the  other  hand,  is  used  when  there  is  no  a  priori  statistical  informa¬ 
tion  is  provided.  With  the  same  AR  model,  both  schemes  are  shown  to  achieve  similar 
performance,  which  is  an  improvement  over  the  Gaussian  assumed  environment. 

D.  Multiuser  Detection  Techniques 

A  narrowband  interference  could  be  a  digital  communication  signal  with  a  data  rate 
much  lower  than  the  spread  spectrum  chip  rate.  This  is  typically  the  case  when  spread 
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spectrum  signals  are  used  in  services  overlaying  on  existing  frequency  band  occupants.  In 
this  case,  the  narrowband  interference  is  a  strong  communication  signal  that  interferes  with 
commercial  DS/SS  communication  systems.  This  type  of  interferers  is  poorly  modeled  as 
either  a  sinusoid  or  an  autoregressive  process.  Because  of  the  similarity  of  the  spread 
spectrum  signal  and  the  digital  interference,  techniques  from  multi-user  detection  theory 
are  applied  to  decode  the  SS  user  signal  and  simultaneously  suppressing  the  interferer  [23]. 

In  order  to  apply  methods  from  multiuser  detection,  the  single  narrowband  interferer  is 
treated  as  a  collection  of  m  spread  spectrum  users,  where  m  is  a  function  of  the  relative 
data  rates  of  the  true  SS  signal  and  the  interference.  That  is,  m  bits  of  the  narrowband 
user  occur  for  each  bit  of  the  SS  user.  As  shown  in  Fig.  4,  and  using  square  waves  for 
illustrations,  each  narrowband  user’s  bit  can  be  thought  as  a  signal  arising  from  a  virtual 
user  having  a  signature  sequence  with  only  one  non-zero  entry.  The  virtual  users  are 
orthogonal,  but  correlated  with  the  SS  user. 

The  optimum  receiver  implementing  the  maximum  likelihood  (ML)  detector  has  a  com¬ 
plexity  that  is  exponential  in  the  number  of  virtual  users,  m.  To  overcome  such  complexity, 
the  optimal  linear  detector  and  decorrelating  detector  are  applied.  While  the  first  requires 
knowledge  of  relative  energies  of  both  the  narrowband  interferer  and  the  SS  user  and  max¬ 
imizes  the  receiver  asymptotic  efficiency,  the  latter  is  independent  of  the  receiver  energies 
and  achieves  the  near-far  resistance  of  the  ML  detector.  The  asymptotic  efficiency  is  the 
limit  of  the  receiver  efficiency  as  the  AWGN  goes  to  zero.  It  characterizes  the  detector 
performance  when  the  dominant  source  of  corruption  is  the  narrowband  interferer  rather 
than  the  AWGN.  The  receiver  efficiency,  on  the  other  hand,  quantifies  the  SS  user  energy 
that  would  achieve  the  same  probability  of  error  in  a  system  with  the  same  AWGN  and 
no  other  users.  The  input  of  both  detectors  is  the  output  of  the  filter  bank  and  consists 
of  filters  matched  to  the  spreading  codes  of  each  active  user,  as  depicted  in  Fig.  5. 

The  following  expressions  are  derived  in  reference  [23]  for  the  probability  of  errors  of  four 
different  detectors.  In  all  cases,  it  is  assumed  that  the  narrowband  signal  is  synchronized 
with  the  SS  signal. 

Conventional  Detector  (CD),  where  the  received  signal  is  sent  directly  to  a  single  filter 
matched  to  the  spreading  code.  The  output  of  the  filter  is  then  compared  to  a  threshold 
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to  yield  the  spread  spectrum  bit  estimate.  This  detector  is  only  optimum  in  the  case  of  a 
single  spread  spectrum  user  in  AWGN.  The  BER  is  given  by 


where  a  —  yjwi/w2,  w\  is  the  received  energy  of  the  narrowband  interference,  w2  is  the 
received  energy  of  the  SS  user  (including  the  process  gain),  p  is  the  vector  formed  by  the 
cross  correlation  between  the  narrowband  interference  waveform  and  the  DS/SS  signal 
waveform,  q  is  the  narrowband  interference  data  bits,  and  {q*}  is  an  ordering  of  the  2m 
possible  values  of  the  vector  of  narrowband  bits. 

Decorrelating  Detector  (DD),  where  the  last  row  of  the  inverse  of  the  cross-correlation 
matrix  of  the  m  +  1  users  is  used  to  multiply  the  output  of  the  m  +  1  matched  filters, 
followed  by  a  threshold  comparison  for  bit  estimate.  The  BER  is  given  by 


Optimum  Linear  Detector  (OLD),  where  the  user  energies  are  used  to  maximize  the 
asymptotic  efficiency.  The  BER  is  given  by 
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where  the  zth  element  of  vector  v  is  given  by 


1  —  Pi  >  Oi 

Vi  —  <- 1  Pi>  a  (29) 

k  pi/  (x  otherwise. 

Ideal  Predictor/Subtracter  (IPS),  which  is  similar  to  the  transversal  filter  excision  tech¬ 
niques  described  in  Section  III-A.  Perfect  knowledge  of  the  narrowband  signal  is  assumed. 
Further,  it  is  assume  that  perfect  prediction  to  the  sample  interior  to  the  narrowband  bit 
is  achieved  and  the  only  error  occurs  when  predicting  at  bit  transitions.  The  expressions 
were  derived  in  [23],  where  one  detector  assumes  zero  bit  estimate  of  the  narrowband  bit 
at  the  transition  and  the  other  detector  takes  this  estimate  to  be  random.  For  the  former 


106 


detector, 


(30) 


where  qp  is  the  estimate  of  q^,  and  p,  defined  only  over  the  chip  interval  emcompassing 
a  narrowband  bit  transition,  is  the  vector  formed  by  the  cross  correlation  between  the 
narrowband  interference  and  the  DS/SS  signal. 

The  performance  of  the  optimum  linear  and  decorrelator  detectors,  representing  the 
multi-user  detection  techniques,  is  shown  in  reference  [23]  to  be  similar  for  different  inter¬ 
ference  power  and  bandwidth.  Both  techniques  significantly  outperform  the  conventional 
detector  and  the  predictor/subtracter,  when  using  a  7  tap  LMS  prediction  filter.  The  im¬ 
provement  in  BER  is  more  pronounced  for  stronger  and  less  narrowband  interferers.  The 
advantage  of  the  decorrelator  detector  over  the  other  conventional  and  adaptive  prediction 
filters  remains  unchanged  when  considering  asynchronous  interference. 

Fig.  6  depicts  the  BER  comparison  between  the  conventional  detector,  decorrelating 
detector,  optimum  linear  detector,  and  ideal  predictor  /subtracter  with  m— 2  and  L  =  63. 
This  figure  is  in  agreement  with  the  performance  figures,  shown  in  [23],  and  conforms  with 
the  same  observations  stated  above.  It  is  clear  from  Fig.  6  that  the  matched  filter  performs 
well  for  weak  interferers.  The  optimum  linear  detector  offers  slight  improvement  over 
the  decorrelating  detector.  The  ideal  predictor /subtracter  outperforms  the  decorrelating 
detector  for  moderate  values  of  interference  power.  It  is  important  to  note,  however,  that 
the  actual  predictor/subtracter  performance  will  have  much  greater  error  probability  [23]. 


E.  Minimum  Mean-Square  Error  Algorithm 

The  minimum  mean-square  error  (MMSE)  algorithm,  originally  proposed  for  suppress¬ 
ing  multiple-access  interference  in  DS/CDMA  multi-user  detection  problems,  has  been 
employed  for  narrowband  interference  mitigation  [24].  Using  the  signal-to-interference  ra¬ 
tio  and  its  upper  bounds  as  a  performance  measure,  the  MMSE  has  been  compared  with 
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linear  and  nonlinear  techniques  in  suppressing  three  types  of  interference,  namely  single¬ 
tone  and  multi-tone  signals,  autoregressive  process,  and  digital  communications  signal 
with  a  data  rate  much  lower  than  the  spread  spectrum  chip  rate.  The  linear  estima¬ 
tors  include  the  conventional  matched  filter  detector,  the  predictor/subtracter,  and  the 
interpolator/subtracter  techniques.  The  latter  is  based  on  using  a  fixed  number  of  past 
and  future  samples  [24].  The  nonlinear  techniques  include  those  based  on  prediction  and 
interpolation.  It  is  shown  that  the  MMSE  detector  completely  suppresses  the  digital  in¬ 
terference,  irrespective  of  its  power,  and  provides  similar  performance  to  the  nonlinear 
interpolator /subtracter  method,  when  dealing  with  AR  type  of  interference. 

IV.  Nonstationary  Interference  Suppression 

The  interference  excision  techniques  discussed  in  the  previous  sections  deal  with  sta¬ 
tionary  or  quasi-stationary  environment.  The  interference  frequency  signature,  or  charac¬ 
teristics,  is  assumed  fixed  or  slowly  time-varying.  None  of  these  techniques  is  capable  of 
effectively  incorporating  the  suddenly  changing  or  evolutionary  rapidly  time-varying  na¬ 
ture  of  the  frequency  characteristics  of  the  interference.  These  techniques  all  suffer  from 
their  lack  of  intelligence  about  interference  behavior  in  the  joint  time-frequency  (t-f)  do¬ 
main  and  therefore  are  limited  in  their  results  and  their  applicability.  For  the  time- varying 
interference  depicted  in  Fig.  7,  frequency-domain  methods  remove  the  frequency  band  A / 
and  ignore  the  fact  that  only  few  frequency  bins  are  contaminated  by  the  interference  at 
a  given  time.  Dually,  time  domain  excision  techniques,  through  gating  or  clipping  the 
interference  over  At,  do  not  account  for  the  cases  where  only  few  time  samples  are  con¬ 
taminated  by  the  interference  for  a  given  frequency.  Applying  either  method  will  indeed 
eliminate  the  interference  but  at  the  cost  of  unnecessarily  reducing  the  desired  signal  en¬ 
ergy.  Adaptive  excision  methods  might  be  able  to  track  and  remove  the  nonstationary 
interference,  but  would  fail  if  the  interference  is  highly  nonlinear  FM  or  linear  FM,  as  in 
Fig.  7,  with  high  sweep  rates.  Further,  the  adaptive  filtering  length  or  block  transform 
length  trades  off  the  temporal  and  the  spectral  resolutions  of  the  interference.  Increasing 
the  step  size  parameter  increases  the  filter  output  errors  at  convergence,  and  causes  an 
unstable  estimate  of  the  interference  waveform. 

The  above  example  clearly  demonstrates  that  nonstationary  interferers,  which  have 
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model  parameters  that  rapidly  change  with  time,  are  particularly  troublesome  due  to  the 
inability  of  single-domain  mitigation  algorithms  to  adequately  ameliorate  their  effects.  In 
this  challenging  situation,  and  others  like  it,  joint  t-f  techniques  can  provide  significant 
performance  gains,  since  the  instantaneous  frequency  (IF),  the  instantaneous  bandwidth, 
and  the  energy  measurement,  in  addition  to  myriad  other  parameters,  are  available.  The 
objective  is  then  to  estimate  the  t-f  signature  of  the  received  data  using  t-f  analysis, 
attenuating  the  received  signal  in  those  t-f  regions  that  contain  strong  interference.  This 
is  depicted  by  the  region  in  between  the  dashed  lines  in  Fig.  7. 

An  FM  interference  in  the  form  u(n)  =  is  solely  characterized  by  its  IF,  which  can 
be  estimated  using  a  variety  of  IF  estimators,  including  the  time-frequency  distributions 
(TFDs)  [25],  [26]. 

The  TFD  of  the  data,  x(n),  at  time  t  and  radian  frequency  u>,  is  given  by 

OO  OO 

Cf(t,u},<f>)  =  l)x(n  +  m  4-  l)x*(x  +  m  —  (32) 

l~ — oo  m~ — oo 

where  )  is  the  time-frequency  kernel  which  is  a  function  of  the  lag  l  and  time-lag 

m.  Several  requirements  have  been  imposed  on  <j>{m ,  l)  to  satisfy  desirable  distribution 
properties,  including  power  localization  at  the  IF.  As  shown  in  eqn.  (32),  the  TFD  is  the 
Fourier  transform  of  a  time-average  estimate  of  the  autocorrelation  function. 

A  time-frequency  notch  filter  can  be  designed,  in  which  the  position  of  the  filter  notch 
is  synchronous  with  the  interference  IF  estimate.  Based  on  the  IF,  two  constraints  should 
exist  to  construct  an  interference  excision  filter  with  desirable  characteristics.  First,  an 
FIR  filter  with  short  impulse  response  must  be  used.  Long  extent  filters  are  likely  to 
span  segments  of  changing  frequency  contents  and,  as  such,  allow  some  of  the  interference 
components  to  escape  to  the  filter  output.  Second,  at  any  given  time,  the  filter  frequency 
response  must  be  close  to  an  ideal  notch  filter  to  be  able  to  null  the  interference  with 
minimum  possible  distortion  of  the  signal.  This  property,  however,  requires  filters  with 
infinite  or  relatively  long  impulse  responses. 

Amin  [27]  has  shown  that  a  linear-phase  five-coefficient  filter  is  effective  in  FM  inter¬ 
ference  excision.  Assuming  exact  IF  values,  the  corresponding  receiver  SINK  is  given 
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by 


SINR  = 


L 

11/8 +  9*2/4' 


(33) 


The  above  expression  shows  that  full  interference  excision  comes  at  the  expense  of  a  change 
in  the  noise  variance  in  addition  of  a  self-noise  form,  as  compared  with  the  non-interference 
case.  The  main  objective  of  any  excision  process  is  to  reduce  both  effects.  The  SINR  in 
(33)  assumes  a  random  IF  with  uniform  distribution  over  [0, 27r].  For  an  interference  with 
fixed  frequency  u>q,  the  receiver  SINR  becomes  dependent  on  loq.  The  receiver  performance 
sensitivity  to  the  interference  frequency  is  discussed  in  detail  in  [27]. 

Wang  and  Amin  [28]  considered  the  performance  analysis  of  the  IF-based  excision  system 
using  a  general  class  of  multiple-zero  FIR  excision  filters  showing  the  dependence  of  the 
BER  on  the  filter  order  and  its  group  delay.  The  effect  of  inaccuracies  in  the  interference 
IF  on  receiver  performance  was  also  considered  as  a  function  of  the  filter  notch  bandwidth. 
Closed  form  approximations  for  SINR  at  the  receiver  are  given  for  the  various  cases. 

One  of  the  drawbacks  to  the  notch  filter  approach  in  [27]  is  the  infinite  notch  depth  due 
to  the  placement  of  the  filter  zeros.  The  effect  is  a  “self-noise”  inflicted  on  the  received 
signal  by  the  action  of  the  filter  on  the  PN  sequence  underlying  the  spread  information 
signal.  This  problem  led  to  the  design  of  an  open-loop  filter  with  adjustable  notch  depth 
based  on  the  interference  energy.  The  notch  depth  is  determined  by  a  variable  embedded 
in  the  filter  coefficients  chosen  as  the  solution  to  an  optimization  problem  which  maximizes 
receiver  SINR.  The  TFD  is  necessary  for  this  work,  even  for  single  component  signals,  be¬ 
cause  simple  IF  estimators  do  not  provide  energy  information.  Amin,  Wang,  and  Lindsey 
accomplished  this  work  in  [29],  incorporating  a  “depth  factor”  into  the  analysis  and  re¬ 
developing  all  the  SINR  calculations.  The  result  was  a  significant  improvement  in  SINR, 
especially  at  mid-range  interference-to-signal  ratios  (ISR’s),  typically  around  0  to  20  dB. 

Instead  of  using  time-varying  excision  filters,  Barbarossa  and  Scaglione  [30]  proposed  a 
two-step  procedure  based  on  dechirping  techniques  commonly  applied  in  radar  algorithms. 
In  the  first  step  the  time  varying  interference  is  converted  to  a  fixed  frequency  sinusoid 
eliminated  by  time  invariant  filters.  The  process  is  reversed.  In  the  second  step  and  the 
interference-free  signal  is  multiplied  by  the  interference  t-f  signature  to  restore  the  DS/SS 
signal  and  noise  characteristics  which  have  been  strongly  impacted  in  the  first  phase. 
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Similar  to  predictor/subtracter  method  discussed  in  Section  III,  Lach,  Amin,  and  Lind¬ 
sey  proposed  synthesis/subtracter  technique  for  FM  interference  using  TFD  [31].  A  replica 
of  the  interference  can  be  synthesized  from  the  t-f  domain  and  subtracted  from  the  incom¬ 
ing  signal  to  produce  an  essentially  interference-free  channel. 

Another  synthesis/subtracter  method  is  introduced  in  [32]  where  the  discrete  evolution¬ 
ary  and  the  Hough  transforms  are  used  to  estimate  the  IF.  The  interference  amplitude  is 
found  by  conventional  methods  such  as  linear  filtering  or  singular  value  decomposition. 
This  excision  technique  applies  equally  well  to  one  or  multi-component  chirp-interferers 
with  constant  or  time-varying  amplitudes  and  with  instantaneous  frequencies  not  neces¬ 
sarily  parametrically  modeled. 

To  overcome  the  drawbacks  of  the  potential  amplitude  and  phase  errors  produced  by 
the  synthesis  methods,  Amin,  Ramineni  and  Lindsey  [33]  proposed  a  projection  filter 
approach  in  which  the  FM  interference  subspace  is  constructed  from  its  t-f  signature. 
Since  the  signal  space  at  the  receiver  is  not  specifically  mandated,  it  can  be  rotated  such 
that  a  single  interferer  becomes  one  of  the  basis  functions.  In  this  way,  the  interference 
subspace  is  one  dimensional  and  its  orthogonal  subspace  is  interference-free.  A  projection 
of  the  received  signal  onto  the  orthogonal  subspace  accomplishes  interference  excision  with 
a  minimal  message  degradation.  The  projection  filtering  methods  compare  favorably  over 
the  previous  notch  filtering  systems. 

In  [34],  Zhang,  Amin,  and  Lindsey  proposed  a  method  to  suppress  more  general  INBI 
signals.  The  interference  subspace  is  constructed  using  t-f  synthesis  methods.  In  different 
to  the  work  in  [31],  the  interferer  is  removed  by  projection  rather  than  subtraction.  To 
estimate  the  interference  waveform,  a  mask  is  constructed  and  applied  such  that  the 
masked  t-f  region  captures  the  interference  energy,  but  leaves  out  most  of  the  DS/SS 
signals. 

Seong  and  Loughlin  have  also  extended  the  projection  method  developed  by  Amin  et. 
al.  [33]  for  excising  constant  amplitude  FM  interferers  from  DS/SS  signals  to  the  case 
of  AM-FM  interferers  [35].  Theoretical  performance  results  (correlator  SNR  and  BER) 
for  the  AM-FM  projector  filter  show  that  FM  estimation  errors  generally  cause  greater 
performance  degradation  than  the  same  level  of  error  in  estimating  the  AM.  The  lower- 
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bound  for  the  correlator  SINR  for  the  AM-FM  projection  filter  for  the  case  of  both  AM 
and  FM  errors  is  given  by 


SINR  = 


l  +  <  +  A 


1  +  a 


2 

A  a 


(* 


where  L  is  the  PN  sequence  length,  A2  is  the  interference  power,  a2  is  the  variance  of 
AWGN,  and  ffia  and  a2A(j)  are  the  variances  of  the  estimation  errors  in  the  AM  and  FM, 
respectively. 

Linear  t-f  signal  analysis  has  also  been  shown  effective  to  characterize  a  large  class  of 
nonstationary  interferers.  Roberts  and  Amin  (36]  proposed  the  use  of  the  discrete  Gabor 
transform  (DGT)  as  a  linear  joint  time-frequency  representation.  The  DGT  can  attenuate 
a  large  class  of  nonstationary  wideband  interferers  whose  spectra  are  localized  in  the  t-f 
domain.  Compared  to  bilinear  TFDs,  the  DGT  does  not  suffer  from  the  crossterm  in¬ 
terference  problems,  and  enjoys  a  low  computational  complexity.  In  [37],  Wei,  Harding, 
and  Bovik  devised  a  DGT-based,  iterative  time-varying  excision  filtering,  in  which  a  hy¬ 
pothesis  testing  approach  was  used  to  design  a  binary  mask  in  the  DGT  domain.  The 
time-frequency  geometric  shape  of  the  mask  is  adapted  to  the  time-varying  spectrum  of 
the  interference.  They  show  that  such  a  statistical  framework  for  the  transform-domain 
mask  design  can  be  extended  to  any  linear  transform.  Both  the  maximum  likelihood  test 
and  the  local  optimal  test  are  presented  to  demonstrate  performance  versus  complexity. 

The  application  of  the  short-time  Fourier  transform  (STFT)  to  nonstationary  interfer¬ 
ence  excision  in  DS/SS  communications  is  considered  in  [38],  [39].  In  both  papers,  due 
to  the  inherent  property  of  STFT  to  trade  off  temporal  and  spectral  resolutions,  several 
STFTs  corresponding  to  different  analysis  windows  were  generated.  In  [38],  Ouyang  and 
Amin  used  a  multiple-pole  data  window  to  obtain  a  large  class  of  recursive  STFTs.  Sub¬ 
sequently,  they  employed  concentration  measures  to  select  the  STFT  that  localizes  the 
interference  in  the  t-f  domain.  This  procedure  is  followed  by  applying  a  binary  excision 
mask  to  remove  the  high-power  t-f  region.  The  remainder  is  synthesized  to  yield  a  DS/SS 
signal  with  improved  signal-to-interference  ratio  (SIR). 

In  [39],  Krongold  et.  al  proposed  multiple  overdetermined  tiling  techniques  and  utilized 
a  collection  of  STFTs  for  the  purpose  of  interference  excision.  Unlike  the  procedure  in  [38], 
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the  authors  in  [39]  removed  the  high-value  coefficients  in  all  generated  STFTs,  and  used 
the  combined  results,  via  efficient  least-square  synthesis,  to  reconstruct  an  interference- 
reduced  signal.  Bultan  and  Akansu  [40]  proposed  a  chirplet-transform-based  exciser  to 
handle  chirp-like  interference  types  in  SS  communications. 

The  block  diagram  in  Fig.  8  depicts  the  various  interference  rejection  techniques  using 
the  time-frequency  methods  cited  above. 

A.  Example 

At  this  point,  in  order  to  further  illustrate  these  excision  methods,  the  work  in  [33]  will  be 
detailed  since  it  includes  comparisons  between  the  two  most  prominent  techniques  based  on 
TFDs  currently  being  studied  —  notch  filtering  and  projection  filtering.  The  signal  model 
is,  as  expected,  given  by  (1),  and  the  major  theme  of  the  work  is  to  annihilate  interference 
via  projection  of  the  received  signal  onto  an  “interference-free”  subspace  generated  from 
the  estimated  interference  characteristics.  This  paper  includes  a  figure,  reprinted  here  as 
Fig.  9,  which  clearly  illustrates  the  trade-offs  between  projection  and  notch  filtering  based 
on  the  ISR.  In  the  legend,  the  variable  a  represents  the  adaptation  parameter  for  the 
notch  filtering  scheme  and  N  represents  the  block  size,  in  samples,  for  a  128  sample  bit 
duration  in  the  projection  method.  Thus,  N  =  128  means  no  block  processing  and  N  =  2 
corresponds  to  64  blocks  per  bit  being  processed  for  projection.  Since  the  projection  and 
non-adaptive  notch  filter  techniques  are  assumed  to  completely  annihilate  the  interference, 
their  performance  is  decoupled  from  the  interference  power,  and  therefore  correctly  indicate 
constant  SINR  across  the  graph.  The  dashed  line  representing  the  notch  filter  with  a  —  0 
is  really  indicating  no  filtering  at  all,  since  the  adaptation  parameter  controls  the  depth 
of  the  notch. 

It  is  evident  from  Fig.  9  that  without  adaptation  a  crossover  point  occurs  around  2  dB 
where  filtering  with  an  infinitely  deep  notch  is  advantageous.  Thus  when  the  interference 
power  exceeds  this  point,  presumably  a  user  would  flip  a  switch  to  turn  on  the  excision 
subsystem.  However,  with  adaptation,  this  process  happens  automatically,  while  giving 
superior  performance  in  the  midrange.  For  the  projection  technique,  the  block  size  de¬ 
termines  receiver  performance  conspicuously  (ceteris  paribus).  Most  important  to  note, 
however,  is  the  superior  performance  of  projection  over  all  methods  when  the  block  size 
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is  equal  to  the  bit  duration,  i.e.  no  block  processing.  It  is  feasible  that  computational 
complexity  may  warrant  a  trade-off  between  SINR  and  block  size,  in  which  case  a  hybrid 
implementation  may  be  of  benefit  —  one  that  automatically  switches  between  adaptive 
notch  filtering  and  projection  depending  on  the  desired  SINR.  In  any  case,  this  example  il¬ 
lustrates  the  parameters  involved  in  the  design  of  modern  excision  filters  for  nonstationary 
interferers. 

V.  Interference  Suppression  for  Frequency  Hopping  Communications 

Interference  rejection  for  FH  is  not  as  well  developed  as  interference  rejection  for  DS  or 
for  CDMA.  In  FH  systems,  the  fast  FH  (FFH)  is  of  most  interest,  and  the  modulation  most 
commonly  used  in  FH  is  frequency-shift  keying  (FSK).  Two  types  of  interference  waveforms 
can  be  categorized,  namely,  partial-band  interference  (PBI)  and  multitone  interference 
(MTI). 

The  effects  of  PBI  and  AWGN  on  several  diversity-combining  receivers  in  FFH/FSK  SS 
communication  systems  have  been  investigated  in  [41],  [42],  [43],  [45].  In  [46],  an  alterna¬ 
tive  method  using  a  fast  Fourier  transform  (FFT)  is  proposed.  In  [41]  the  authors  present 
an  automatic  gain-control  (AGC)  receiver  using  a  diversity  technique.  In  this  method, 
each  soft-decision  square-law  detected  MARK  and  SPACE  filter  output  is  weighted  by 
the  inverse  of  the  noise  power  in  the  slot  prior  to  linear  combining.  This  method  is  near- 
optimal  (in  terms  of  SNR)  if  the  exact  information  of  noise  and  interference  power  can 
be  obtained.  A  similar  clipped-linear  combining  receiver  was  also  reported  in  [43].  Due 
to  the  difficulty  of  such  information,  self-normalizing  receivers  [42]  and  the  ratio-statistic 
combining  technique  [44]  use  the  output  values  of  the  square-law  detector  in  each  hop  to 
derive  a  weight  or  normalizing  factor.  The  performance  of  these  two  methods  is  shown  to 
be  comparable  to  that  of  the  square-law  clipper  receiver. 

In  [45]  the  authors  describe  an  FFH  receiver  that  employs  a  prewhitening  filter  to  reject 
NBI.  For  binary  FSK  modulations,  it  is  shown  that  the  FFH  signal  is  statistically  uncorre¬ 
lated  at  lag  values  of  T/,/(4Ar/l)  where  Th  is  the  hop  duration  and  2Nh  is  the  total  number 
of  frequency  slots  (i.e.,  there  are  Nh  MARK  and  Nh  SPACE’s).  Thus,  as  in  the  DS  case, 
NBI  can  be  predicted  and  suppressed  independently  of  the  desired  FFH  signal.  Using  the 
complex  LMS  algorithm  to  update  the  prewhitening  filter  coefficients,  this  technique  is 
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shown  to  compare  favorably  with  the  maximal-ratio  combiner  diversity  technique.  When 
the  interferer  is  wide-sense  stationary,  the  prewhitening  filter  based  receiver  provides  per¬ 
formance  approaching  that  of  the  AGC  receiver  and  at  least  2-3  dB  superior  to  that  of 
the  self-normalizing  receiver.  However,  when  hostile  interference  is  present,  the  adaptive 
prewhitening  filter  technique  may  not  be  able  to  track  the  interference  rapidly  enough.  In 
this  case,  nonparametric  techniques  such  as  the  self-normalizing  receiver  must  be  used  to 
reject  the  jammed  hops. 

Reed  and  Agee  [47]  have  extended  and  improved  on  the  idea  of  whitening  by  using  a 
time-dependent  filter  structure  to  estimate  and  remove  interference,  based  on  the  interfer¬ 
ence  spectral  correlation  properties.  In  this  method,  the  detection  of  FH/SS  in  the  presence 
of  spectrally  correlated  interference  is  nearly  independent  of  the  SIR.  The  process  can  be 
viewed  as  a  time-dependent  whitening  process  with  suppression  of  signals  that  exhibit  a 
particular  spectral  correlation.  The  technique  is  developed  from  the  maximum-likelihood 
estimate  of  the  spectral  frequency  of  a  frequency  agile  signal  received  in  complex  Gaussian 
interference  with  unknown  spectral  correlation.  The  resulting  algorithm  uses  the  corre¬ 
lation  between  spectrally  separated  interference  components  to  reduce  the  interference 
content  in  each  spectral  bin  prior  to  the  whitening/detection  operation. 

In  [46],  an  alternative  approach  to  suppress  PBI  using  the  FFT  is  proposed.  The  major 
attraction  of  FFT-based  implementation  lies  in  the  ability  to  achieve  guaranteed  accuracy 
and  perfect  reproducibility. 

For  suppression  of  MTI,  basically  the  same  processing  methods  applied  for  PBI  can  be 
employed.  However,  the  performance  analyses  differ  from  those  for  PBI  situations.  The 
performance  depends  on  the  distribution  of  the  MTI  and,  in  turn,  how  many  bands  are 
contaminated  by  MTI.  Performance  analyses  of  FFH  SS  systems  are  presented  in  [48],  [49] 
for  linear  combining  diversity,  in  [50]  clipped  diversity,  in  [51],  [52]  for  maximum  likelihood 
and  product-combining  receivers. 
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Fig.  3  Trans form- domain  notch  filtering. 
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Fig.  4  Virtual  CDMA  systems  (synchronous  case) . 
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Multiuser  detector  structure. 


Excision  methods  for  nonstationary  interferers. 
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Fig. 7  Interference  rejection  techniques 


Fig.  8  Comparison  between  projection  and  notch  filtering  excision 
methods . 


123 


In  B.  Boashash,  Time-Frequency  Signal  Analysis  and  Processing: 

A  Comprehensive  Reference ,  Elsevier,  Oxford,  UK,  2003. 

Spatial  Time-Frequency  Distributions 
and  Their  Applications 

Moeness  G.  Amin  and  Yimin  Zhang 

Department  of  Electrical  and  Computer  Engineering,  Villanova 
University,  Villanova,  PA  19085,  USA 

I.  Spatial  Time-Frequency  Distributions 

The  evaluation  of  quadratic  time-frequency  distributions  (TFDs)  of  nonstationary  sig¬ 
nals  impinging  on  a  multi-sensor  receiver  yields  spatial  time-frequency  distributions  (STFDs), 
which  permit  the  application  of  eigenstructure  subspace  techniques  to  solving  a  large  class 
of  channel  estimation  and  equalization,  blind  source  separation  (BSS),  and  high  resolution 
direction-of-arrival  (DOA)  estimation  problems  [1],  [2],  [3].  STFD  based  techniques  are 
appropriate  to  handle  sources  of  nonstationary  waveforms  that  are  highly  localized  in  the 
time-frequency  (t-f)  domain.  In  the  area  of  BSS,  the  use  of  the  STFDs  allows  the  separa¬ 
tion  of  sources  with  identical  spectral  shape,  but  with  different  t-f  localization  properties, 
i.e.,  different  t-f  signatures.  For  both  source  separation  and  DOA  estimation  problems, 
spreading  the  noise  power  while  localizing  the  source  energy  in  the  t-f  domain  amounts  to 
increasing  the  robustness  of  eigenstructure  signal  and  noise  subspace  estimation  methods 
with  respect  to  channel  and  receiver  noise.  This  in  turn  leads  to  an  improvement  of  spatial 
resolution  and  source  separation  performance. 

The  quadratic  class  of  STFD  matrix  of  a  signal  vector  x(f)  is  defined  as 

D  xx(tJ)=  f  f  f  g{v,T)x(u  +  ^)xH(u-^)ejMvu-vt-fT)drdudv,  (1) 

J  —  00  J  —  OO  J — oo  Z  Z 

where  g(v,r )  is  the  kernel  function. 

In  narrowband  array  processing,  when  n  signals  arrive  at  an  m-element  array  (see  Fig. 

I),  the  linear  data  model 

x(f)  =  y  (t)  +  n(t)  =  Ad  (t)  +  n  (t)  (2) 
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is  commonly  assumed,  where  x(t)  is  the  mxl  data  vector  received  at  the  array,  d(t)  is  the 
n  x  1  source  data  vector,  the  m  x  n  spatial  matrix  A  =  [ai  •  •  •  an]  represents  the  mixing 
matrix,  a*  is  the  steering  vector  of  ith  signal,  and  n(t)  is  an  additive  noise  vector  whose 
elements  are  modeled  as  stationary,  spatially  and  temporally  white,  zero-mean  complex 
random  processes,  independent  of  the  source  signals. 


ra-element  array  with  n  signal  arrivals. 


Under  the  uncorrelated  signal  and  noise  assumption  and  the  zero-mean  noise  property, 
the  expectation  of  the  crossterm  TFD  matrices  between  the  signal  and  noise  vectors  is 
zero,  i.e.,  E  [Dyn(t,  /)]  =  E  [Dny(t,  /)]  =  0,  and  it  follows 

E  [Dxx(i,  /)]  =  Dyy  (t,  f )  +  E  [Dnn(f,  /)]  =  ADdd(i,  f)AH  +  a2 1,  (3) 

where  a2  is  the  noise  power,  and  I  is  the  identity  matrix.  Equation  (3)  is  similar  to 
that  which  has  been  commonly  used  in  array  processing  based  on  second-order  statis¬ 
tics,  relating  the  signal  correlation  matrix  to  the  data  spatial  correlation  matrix  [1],  This 
implies  that  key  problems  in  various  applications  of  array  processing,  specifically  those 
dealing  with  nonstationary  signal  environments,  can  be  approached  using  quadratic  trans¬ 
formations.  If  Ddd(t, /)  is  a  full-rank  matrix,  the  two  subspaces  spanned  by  the  principal 
eigenvectors  of  D xx(t,f)  and  the  columns  of  A  become  identical.  In  this  case,  direction 
finding  techniques  based  on  eigenstructures  can  be  applied.  If  Ddd(i,  /)  is  diagonal,  i.e., 
the  signal  cross-TFDs  at  the  t-f  point  ( t ,  /)  are  zeros,  then  both  the  mixing  matrix  and 
the  signal  waveforms  can  be  recovered  using  BSS  methods. 
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II.  Fundamental  Properties 

There  are  five  key  advantages  of  array  signal  processing  using  STFD.  In  order  to  properly 
explain  these  advantages,  we  use  the  diagram  in  Fig.  II.  We  consider  two  sources  A  and 
B  incident  on  a  multi-sensor  array.  Source  A  occupies  the  t-f  region  Ra,  whereas  source  B 
occupies  the  t-f  region  Rb .  The  t-f  signatures  of  the  two  sources  overlap,  but  each  source 
still  has  a  t-f  region  that  is  not  intruded  over  by  the  other  source. 


Signals  with  different  time-frequency  signature. 

1)  Equation  (3)  can  be  easily  derived  for  any  arbitrary  joint- variables.  Time  and  fre¬ 
quency  are  indeed  the  two  most  commonly  used  and  physically  understood  parameters. 
However,  by  replacing  the  STFDs  by  spatial  arbitrary  joint- variable  distributions,  one 
can  relate  the  sensor  joint-variable  distributions  to  the  sources  joint-variable  distributions 
through  the  same  mixing  matrix  A.  As  shown  in  the  Examples  Section,  there  are  situa¬ 
tions  where  it  is  preferable  to  consider  other  domains  such  as  the  ambiguity  lag-Doppler 
domain,  where  the  locations  of  the  signals  and  their  cross-terms  are  guided  by  properties 
and  mechanisms  different  than  those  associated  with  the  t-f  domain. 

2)  Equation  (3)  is  valid  for  all  t-f  points.  It  is  well  known  that  direction  finding  tech¬ 
niques  require  Ddd(£,  /)  to  be  full  rank,  preferably  diagonal.  On  the  other  hand,  BSS 
techniques  demand  the  diagonal  structure  of  the  same  matrix  without  degenerate  eigen¬ 
values.  These  properties  along  with  high  signal-to-noise  ratio  (SNR)  requirements  may  be 
difficult  to  achieve  using  a  single  t-f  point.  Two  different  methods  can  be  used  for  inte¬ 
grating  several  t-f  points  into  equation  (3).  One  method  is  based  on  a  simple  averaging 
performed  over  the  signatures  of  the  sources  of  interest,  whereas  the  second  method  is 
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based  on  incorporating  several  desired  t-f  points  into  joint  diagonalization  or  joint  block- 
diagonalization  schemes. 

3)  The  TFD  of  the  white  noise  is  distributed  all  over  the  t-f  domain,  whereas  the  TFDs 
of  the  source  waveforms  are  likely  to  be  confined  to  much  smaller  regions.  Referring  to 
Fig.  II,  the  noise  is  spread  over  both  Ra  and  Rb  as  well  as  the  complement  region  Rc. 
If  the  t-f  points  ( t ,  /)  used  in  either  the  averaging  or  joint  diagonalization  approaches 
belong  to  the  noise  only  region  Rc,  then  no  information  of  the  incident  waveforms  is  used 
and,  as  such,  no  reasonable  source  localization  and  signal  separation  outcomes  can  be 
obtained.  On  the  other  hand,  if  all  points  ( t ,  /)  in  Fig.  II  are  used,  and  the  employed  TFD 
satisfies  the  marginal  constraints,  then  it  can  be  easily  shown  that  only  the  signal  average 
power  is  considered.  As  a  result,  the  problem  simplifies  to  the  second-order  covariance 
based  matrix  approach,  traditionally  used  in  high  resolution  DOA  estimation.  This  is 
an  important  feature,  as  it  casts  the  conventional  techniques  as  special  cases  of  the  array 
signal  processing  framework  based  on  t-f  analysis.  Finally,  if  we  confine  the  ( t ,  /)  points  to 
Ra  and  Rb,  then  only  the  noise  part  in  these  regions  is  included.  The  result  of  leaving  out 
the  points  ( t ,  /)  that  are  not  part  of  the  t-f  signatures  of  the  signal  arrivals  is  enhancing  the 
input  SNR,  which  is  utilized  by  the  source  localization  and  signal  separation  techniques. 

4)  By  only  selecting  t-f  points  that  belong  to  the  t-f  signature  of  one  source,  then 
this  source  will  be  the  only  one  considered  by  equation  (3).  This  selection,  in  essence, 
is  equivalent  to  implicitly  performing  spatial  filtering  and  removing  other  sources  from 
consideration.  It  is  important  to  note,  however,  that  such  removal  does  not  come  at 
the  expense  of  reduction  of  the  number  of  degrees-of-freedom  (DOFs),  as  it  is  the  case  in 
beamspace  processing,  but  the  problem  remains  a  sensor  space  processing  with  the  original 
number  of  DOFs  kept  intact.  This  property  represents  a  key  contribution  of  TFDs  to  the 
direction  finding  and  DOA  estimation  areas.  An  antenna  array  can  be  used  to  localize 
a  number  of  sources  equal  or  even  greater  than  its  number  of  sensors.  The  fundamental 
condition  is  that  there  must  be  t-f  regions  over  which  the  respective  t-f  signatures  of  the 
sources  do  not  overlap.  Referring  to  Fig.  II  and  considering  the  case  of  two  sensors,  if 
all  t-f  points  incorporated  in  direction  finding  belong  to  region  Ra  and  not  Rb,  then  the 
signal  subspace  defined  by  equation  (3)  is  one-dimensional.  Thus,  by  excluding  source 
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B,  a  one-dimensional  noise  subspace  is  established.  This  allows  us  to  proceed  with  high 
resolution  techniques  for  localization  of  source  A.  In  a  general  scenario,  one  can  localize 
one  source  at  a  time  or  a  set  of  selected  sources,  depending  on  the  array  size,  overlapping 
and  distinct  t-f  regions,  and  the  dimension  of  the  noise  subspace  necessary  to  achieve  the 
required  resolution  performance.  The  same  concepts  and  advantages  of  t-f  point  selection 
discussed  above  for  direction  finding  can  be  applied  to  BSS  problems. 

5)  The  a  priori  knowledge  of  some  temporal  characteristics  or  the  nature  of  time- varying 
frequency  contents  of  the  sources  of  interest  may  permit  us  to  directly  select  the  t-f  regions 
used  in  equation  (3).  For  instance,  it  is  known  that,  in  the  ambiguity  domain,  all  fixed 
frequency  sinusoidal  signals  map  to  the  time-lag  axis.  By  only  incorporating  the  points 
on  this  axis,  we  have,  in  fact,  opted  to  separate  and  localize  all  narrowband  signals  in 
broadband  communications  platforms. 

III.  Examples 

In  this  Section,  we  present  simulation  examples  to  demonstrate  the  fundamental  offer¬ 
ings  discussed  in  the  previous  Section.  Time-frequency  MUSIC  (t-f  MUSIC),  ambiguity- 
domain  MUSIC  (AD-MUSIC),  and  the  BSS  based  on  STFDs  are  three  different  techniques 
chosen  for  the  demonstration.  The  algorithms  involved  in  the  implementation  of  the  tech¬ 
niques  are  given  in  Tables  I,  II,  and  III  [2],  [4],  [1]. 

Example  I  [4].  Consider  the  scenario  of  a  four-element  equi-spaced  linear  array  spaced 
by  half  a  wavelength,  where  one  chirp  signal  and  two  sinusoidal  signals  are  received.  The 
data  record  has  128  samples.  All  three  signals  have  the  same  SNR  of  20  dB.  The  DOAs 
of  the  chirp  signal  and  the  two  sinusoidal  signals  are  15°,  10°,  and  0°,  respectively.  While 
the  ambiguity  function  of  the  chirp  signal  sweeps  the  ambiguity  domain  with  contribution 
at  the  origin,  the  exact  autoterm  ambiguity  function  of  the  narrowband  arrivals  Si(t)  and 
s2(t)  is  zero  for  non-zero  frequency-lags  and  may  have  non-zero  values  only  along  the 
vertical  axis  v  =  0. 

In  this  simulation  example,  we  selected  24  points  on  the  time-lag  axis,  excluding  the 
origin,  and  as  such  emphasizing  the  narrowband  components.  Fig.  Ill  shows  the  ambiguity 
function  where  the  two  vertical  lines  away  from  the  origin  represent  the  crossterms  between 
the  sinusoidal  components.  Fig.  Ill  shows  the  two  estimated  spatial  spectra  for  three 
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Time- Frequency  MUSIC. 


STEP  I  Form  K  matrices  Dxx(t;,  ft)  for  the  selected  (ti,  fo)  points, 
i— 

STEP  II  The  eigenvectors  of  E  [Dxx(t,  /)]  corresponding  to  the  m  -  n 
smallest  eigenvalues, ei,  •  •  • ,  em_n,  are  obtained  by  joint  block- 
diagonalization,  or  the  eigen-decomposition  of  averaged  matrix 

ilDJfe/,). 

^  i=l 

STEP  III  Estimate  the  number  of  signals  from  the  eigenvalues,  and 

estimate  the  DOAs  from  the  peaks  of  the  t-f  MUSIC  spectra 
f(6)  =  |E^a(0)|  2,  where  E„  =  [ei,  •  •  •  ,em_n],  and  a(0 )  is 
the  steering  vector  corresponding  to  DOA  6. 

Ambiguity-Domain  MUSIC 


Ambiguity-Domain  MUSIC  follows  the  same  procedure  as  time-frequency 
MUSIC  by  using  Dxx(uj,  7*)  instead  of  Dxx(tj,  /*),  i  =  1,  •  •  • ,  K. 


independent  trials,  one  corresponds  to  the  conventional  method  and  the  other  corresponds 
to  the  AD-MUSIC.  There  are  two  dominant  eigenvalues  for  the  case  of  the  AD-MUSIC, 
since  we  have  not  deliberately  considered  the  chirp  signal  through  our  careful  selection  of 
the  ambiguity-domain  points.  It  is  clear  that  the  AD-MUSIC  resolves  the  two  sinusoidal 
signals,  while  the  conventional  MUSIC  could  not  separate  the  three  signals. 

Example  II  [5].  Consider  a  uniform  linear  array  of  eight  sensors  separated  by  half  a  wave¬ 
length.  Two  chirp  signals  emitted  from  two  sources  positioned  at  (di,02)  =  (—10°,  10°), 


129 


Blind  Source  Separation  Based  on  STFDs 


STEP  I 

Estimate  the  auto-correlation  matrix  Rxx  from  T  data  samples. 

Denote  by  Ai,  •  •  • ,  Xn  the  n  largest  eigenvalues  and  hi,  •  •  • ,  h„  the 

corresponding  eigenvectors  of  Rxx. 

STEP  II 

An  estimate  a2  of  the  noise  variance  is  the  average  of  the  m  —  n 

smallest  eigenvalues  of  RxX.  The  whitening  matrix  is  formed  as 

w  =  [(Ai  -  d2)-5h1;  • .  • ,  (A„  -  <72)-!hn]H. 

STEP  III 

Form  K  matrices  by  computing  the  STFD  of  whitened  vector 

z (t)  =  Wx(f)  for  a  fixed  set  of  (t;,  /*)  points,  *  =  1,  •  •  • ,  K, 

corresponding  to  signal  autoterms. 

STEP  IV 

A  unitary  matrix  U  is  then  obtained  as  joint  diagonalizer  of  the 

set  Dzz(ij,  fi ),  i  —  1,  *  *  • ,  K . 

STEP  V 

The  source  signals  are  estimated  as  s (t)  —  UHWx(f),  and  the 
mixing  matrix  A  is  estimated  as  A  =  W#U. 

respectively.  The  data  record  has  1024  samples.  The  start  and  end  frequencies  of  the 
chirp  signal  of  the  source  at  9i  are  fsi  =  0  and  fe i  —  0.5,  while  the  corresponding  two 
frequencies  for  the  signal  of  the  other  source  at  92  are  fs2  =  0.5  and  fe2  =  0,  respectively. 

Fig.  Ill  displays  the  standard  deviations  of  the  DOA  estimation  9\  versus  SNR.  The 
curves  in  this  figure  show  the  theoretical  and  experimental  results  of  the  conventional 
MUSIC  and  t-f  MUSIC.  Pseudo  Wigner-Ville  distribution  with  window  length  L  =  33 
and  129  are  considered.  The  Cramer-Rao  Bound  (CRB)  is  also  shown  in  Fig.  III.  Both 
signals  are  selected  when  performing  t-f  MUSIC.  Simulation  results  are  averaged  over  100 
independent  trials  of  Monte  Carlo  experiments.  The  advantages  of  t-f  MUSIC  in  low  SNR 
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V 


The  ambiguity  functions  of  the  chirp  signal  and  two  sinusoidal  signals. 


AD-MUSIC 


The  estimated  spatial  spectra  of  AD-MUSIC  and  conventional  MUSIC. 


cases  are  evident  from  this  figure.  Fig.  Ill  shows  estimated  spatial  spectra  at  SNR=— 20 
dB  based  on  t-f  MUSIC  ( L  =  129)  and  the  conventional  MUSIC.  The  t-f  MUSIC  spectral 
peaks  are  clearly  resolved. 

Example  III  [1].  In  Fig.  Ill,  we  show  an  example  of  the  application  of  STFDs  to  the  BSS 
problem.  A  three-element  equi-spaced  linear  array  is  considered  where  the  interelement 
spacing  is  half  a  wavelength.  Two  chirp  signals  arrive  at  —10°  and  10°,  respectively.  The 
number  of  data  samples  used  to  compute  the  STFD  is  128.  The  number  of  t-f  points 
employed  in  the  joint  diagonalization  is  p=128,  with  equal  number  of  points  on  each 
signature.  Fig.  111(b)  shows  the  Choi-Williams  distributions  of  two  linear  mixtures  of  the 
original  chirp  signals  depicted  in  Fig.  111(a),  corresponding  to  the  data  at  the  first  and  the 
second  sensors.  Using  the  STFDs,  we  are  able  to  recover  the  original  signals  from  their 
observed  mixture,  as  shown  in  Fig.  III(c). 
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IV.  Crossterm  Issues  in  STFD 


There  are  two  sources  of  crossterms.  The  first  type  are  the  crossterms  that  are  the  results 
of  the  interactions  between  the  components  of  the  same  source  signal.  The  other  type 
of  crossterms  are  those  generated  from  the  interactions  between  two  signal  components 
belonging  to  two  different  sources.  These  crossterms  are  associated  with  cross-TFDs  of  the 
source  signals  and,  at  any  given  t-f  point,  they  constitute  the  off-diagonal  entries  of  the 
source  TFD  matrices  Dja  (t, /)  defined  in  (3).  Although  the  off-diagonal  elements  do  not 
necessarily  destroy  the  full-rank  matrix  property  necessary  for  direction  finding  application 
[6],  they  violate  the  basic  assumption  in  the  problem  of  source  separation  regarding  the 
diagonal  structure  of  the  source  TFD  matrix.  We  must  therefore  select  the  t-f  points  that 
belong  to  autoterm  regions  where  crossterm  contributions  are  at  minimum,  e.g.,  by  using 
a  priori  information  of  the  source  signals. 

The  method  of  spatial  averaging  of  the  STFD  introduced  in  [7]  does  not  reduce  the 
crossterms  as  in  the  case  with  reduced  interference  distribution  kernels,  but  rather  move 
them  from  their  locations  on  the  off-diagonal  matrix  entries  to  be  part  of  the  matrix 
diagonal  elements.  The  other  parts  of  the  matrix  diagonal  elements  represent  the  contri¬ 
bution  of  the  autoterms  at  the  same  point.  Therefore,  not  only  we  are  able  to  set  the 
off-diagonal  elements  of  the  source  TFD  matrix  to  zeros,  but  also  we  can  improve  per¬ 
formance  by  selecting  the  t-f  points  of  peak  values,  irrespective  of  whether  these  points 
belong  to  autoterm  or  crossterm  regions. 


The  standard  deviations  of  DOA  estimation  6\  vs.  SNR. 
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A.  Summary  and  Conclusions 

The  spatial  time-frequency  distribution  (STFD)  is  an  important  tool  for  temporal  and 
spatial  separations  of  sources  emitting  nonstationary  signals.  It  is  a  discriminatory  tool 
that  allows  a  consideration  of  only  a  subset  of  source  signals  impinging  on  a  multi-sensor 
receiver.  This  property  enhances  signal  parameter  estimation  and  permits  direction  finding 
and  signal  separation  to  be  applied  to  a  number  of  sources  that  is  equal  or  even  exceeds 
the  number  of  sensors. 

All  material  presented  in  this  article  is  based  on  the  model  (2).  One  important  change  in 
the  direction  of  the  research  in  the  time-frequency  array  signal  processing  area  was  given  in 
[8],  where  the  strict  model  of  (2)  was  relaxed  and  a  direction  finding  technique  employing 
a  STFD-based  wideband  root-MUSIC  was  proposed.  Another  research  direction  is  the 
utilization  and  integration  of  crossterms  into  STFDs.  It  has  recently  been  shown  [9]  that 
source  separation  can  be  performed  based  on  both  autoterms  and  crossterms  through  joint 
diagonalization  and  joint  anti-diagonalization  schemes  of  STFD  matrices. 
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(a)  TFDs  of  the  source  signals 
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(b)  TFDs  of  the  mixed  signals 


TFO  of  separated  signal : 


(c)  TFDs  of  the  separated  signals 
Blind  source  separation  based  on  STFDs. 
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Abstract 

Time-frequency  distributions  (TFDs)  have  evolved  to  be  a  powerful  technique  for  nonstationaiy 
signal  analysis  and  synthesis.  With  the  use  of  a  multi-sensor  array,  spatial  time-frequency 
distributions  (STFDs)  have  been  developed  and  successfully  applied  to  high-resolution  direction-of- 
arrival  (DO A)  estimations  and  blind  recovery  of  the  source  waveforms.  In  this  paper,  we  introduce  the 
spatial  polarimetric  time-frequency  distribution  (SPTFD)  as  a  platform  to  process  nonstationary  array 
signals  with  two  orthogonal  polarization  components,  such  as  horizontal  and  vertical.  The  use  of  dual 
polarization  empowers  the  STFDs  with  additional  degrees-of-freedom  (DOFs)  and  improves  the 
robustness  of  the  signal  and  noise  subspaces.  This  improvement  serves  to  enhance  DOA  estimation 
and  signal  recovery.  To  demonstrate  the  effectiveness  of  the  SPTFD  platform,  the  polarimetric  time- 
frequency  ESPRIT  (PTF-ESPRIT)  method  is  proposed  and  is  shown  to  outperform  time-frequency, 
polarimetric,  and  conventional  ESPRIT  methods. 


1  This  work  was  supported  in  part  by  the  ONR  under  Grant  No.  N00014-98-1-0176  and  DARPA  under  Grant 
No.  MDA972-02- 1-0022.  The  content  of  the  information  does  not  necessarily  reflect  the  position  or  policy  of 
the  Government,  and  no  official  endorsement  should  be  inferred. 


136 


1  Introduction 


Over  the  past  two  decades,  time-frequency  distributions  (TFDs)  have  evolved  to  be  a  powerful 
technique  for  nonstationary  signal  analysis  and  synthesis  in  the  areas  of  speech,  biomedicine, 
automotive  industry,  and  machine  monitoring  [1-5].  In  radar  signal  processing,  the  time-frequency 
signal  representation,  in  its  linear  and  bilinear  forms,  has  been  used  in  target  detection,  classification, 
and  clutter  suppression  [6-10].  Most  recently,  the  spatial  dimension  has  been  incorporated,  along  with 
the  time  and  frequency  variables,  into  quadratic  and  higher-order  time-frequency  distributions,  and 
led  to  the  introduction  of  spatial  time-frequency  distributions  (STFDs)  for  sensor  signal  processing. 
The  STFD  has  been  successfully  applied  to  high-resolution  direction-of-arrival  (DOA)  estimations 
and  blind  recovery  of  the  source  waveforms  impinging  on  a  multi-sensor  receiver,  specifically  those 
of  nonstationary  temporal  characteristics  [11-19]. 

Polarization  and  polarization  diversity,  on  the  other  hand,  have  been  proven  to  be  very  effective  in 
wireless  communications  and  various  types  of  radar  systems  [20,21].  Antenna  and  target  polarization 
properties  are  widely  employed  in  remote  sensing  and  synthesis  aperture  radar  (SAR)  applications 
[22-24].  Presently,  airborne  and  spacebome  platforms  as  well  as  meteorological  radars  include 
polarization  information  [25,26].  Additionally,  polarization  plays  an  effective  role  for  target 
identification  in  the  presence  of  clutter  [27,28].  Polarization  has  also  been  incorporated  in  array 
antennas  for  improved  signal  parameter  estimation,  such  as  DOA  and  time-of-arrival  (TOA)  [29-31]. 

Despite  the  extensive  research  work  performed  in  time-frequency  signal  representations  and 
polarimetric  signal  processing  methods,  these  two  important  areas  have  not  been  coupled  or 
considered  within  the  same  platform.  In  numerous  communications  and  radar  applications,  moving 
sources/targets  often  generate  time-varying  Doppler  frequency,  with  well  defined  Doppler-frequency 
signatures.  Also,  mechanical  vibration  or  rotation  of  structure  in  a  source/target  typically  includes 
frequency  modulation  on  returned  signals.  Time-frequency  methods  have  been  proposed  to 
characterize  the  Doppler-frequency  signature  as  well  as  to  analyze  micro-Doppler  phenomenon. 
However,  little  attention,  if  any,  has  been  paid  to  the  fact  that  the  return  signal  from  a  moving  or 
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vibrating  target  may  have  a  time-varying  polarization.  The  polarimetric  Doppler  frequency  signature 
contains  valuable  information  that  the  single-polarization  Doppler  frequency  signature  fails  to 
provide. 

In  this  paper,  we  develop  the  spatial  polarimetric  time-frequency  technique  for  multi-antenna 
receivers.  This  technique  utilizes  not  only  the  time-varying  Doppler  frequency  signatures,  but  also  the 
polarization  signatures,  whether  they  are  stationary  or  time-varying.  The  signal  polarization 
information  empowers  the  spatial  time-frequency  distributions  (STFDs),  as  it  retains  the  integrity  of 
eigenstructure  methods  and  improves  the  robustness  of  the  respective  signal  and  noise  subspaces 
under  low  signal-to-noise  ratio  (SNR)  and  coherent  signal  environment. 

The  importance  of  this  technique  stems  from  the  fact  that  targets,  emitters,  and  scatters,  when 
changing  their  positions,  are  likely  to  produce  Doppler  and  polarization  profiles  that  are  time- 
dependent.  If  the  field  of  view  covers  multiple  sources,  then  target  detection,  source  location,  and 
signal  recovery  benefit  from  distinctions  in  polarization,  spatial,  and  time-frequency  signatures.  With 
polarization  no  longer  decoupled  from  the  signal  time-varying  spectrum  characteristics,  high 
resolution  imaging  and  DOA  estimation  can  be  achieved  over  the  cases  where  the  decoupling  is 
enforced. 

The  focus  of  this  paper  is  limited  to  the  proposition  of  the  SPTFD  and,  as  its  application  example  to 
demonstrate  the  effectiveness  of  this  technique,  polarimetric  time-frequency  ESPRIT  (PTF-ESPRIT) 
for  DOA  estimation  of  noncoherent  sources  is  considered  [32].  The  application  to  a  MUSIC-like 
method  for  non-coherent  and  coherent  signals  is  introduced  separately  in  [33-35]. 

This  paper  is  organized  as  follows.  Section  2  introduces  the  signal  model  and  briefly  reviews  the 
definition  of  STFD.  In  Section  3  the  polarimetric  time-frequency  distribution  (PTFD)  and  SPTFD  are 
introduced  and  defined.  The  PTF-ESPRIT  is  proposed  in  Section  4.  Section  5  presents  simulation 
results  and  Section  6  concludes  this  paper. 
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Table  1.  Notations 


A  Matrix 

a  Vector 

(.)r  Transpose 


(.)*  complex  conjugate 

(.)[,]  polarization  index 

(.)tf  Hermitian 


2  Signal  Model 

2.1  Time-Frequency  Distributions 

The  Cohen's  class  of  TFDs  of  a  signal  *(0  is  expressed  as  [4] 

/)„(/,/)=  -u,r)x(t  +  r/2)x\t -r/2)  e~i2*fTdudr  0) 

where  t  and /represent  the  time  index  and  frequency  index,  respectively.  The  kernel  (pit,  t)  uniquely 
defines  the  TFD  and  is  a  function  of  the  time  and  lag  variables.  In  this  paper,  all  the  integrals  are  from 

—  oo  to  60  . 

The  TFDs,  given  in  Eq.  (1),  have  been  shown  to  be  a  powerful  tool  in  analysis  of  signals  with  time- 
varying  frequencies.  They  are  used  in  different  applications  for  ninstationary  signal  detection, 
classification,  and  discriminations  [1-5].  The  cross-sensor  TFD  of  two  signals  Xi(r)  and  x2(t )  is  defined 
by 

Dv  ( t ,  f)  =  J"jV(r  -  u, t) x,  (t  +  r/2 )x2(t  - r/2)  e'^^dudr .  (2) 


2.2  Spatial  Time-Frequency  Distributions 

Consider  a  narrowband  array  processing  problem,  where  n  signals  arrive  at  an  m-element  array. 
The  following  linear  data  model  is  assumed, 
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x(r)  =  As(0  +  n  (0 


(3) 


where  the  mxn  spatial  matrix  A  =  [a,,a2,—,aJ  is  the  mixing  matrix  which  holds  the  spatial 
information.  In  this  paper,  structured  mixing  matrix,  that  is,  A  =  A(0)  =  [a(#,),a(02),  ••■,a(0n)]  is 
assumed,  where  ©=[<?,,  02,  •••,#„]  .  The  elements  of  the  mx  1  vector  x(f)  ,  which  represents  the 
measured  or  sensor  data,  are  multi-component  signals,  while  the  elements  of  the  nxl  vector 
s(f )  =  [s, (f ).  s2  (/),-•-,  s„ (/)]r  are  often  mono-component  signals.  n(t)  is  an  mxl  additive  noise  vector, 
which  consists  of  independent  zero-mean,  white  and  Gaussian  distributed  processes. 

The  SIT'D  of  a  data  vector  x(t)  is  expressed  as  [1 1] 

D„(r,/)  =  jj^(r  -  u,  r)x(r  +  t/2)xh  (t  -  t/2)  e~n*ftdudr  (4) 

where  the  (i,A:)th element  of  D^O,/)  is  defined  as  [Dxx0,/)1*  =  DVt{t,f  ) ,  i,k  =1,  2,  •••,  m. 

Due  to  the  linear  data  model,  the  noise-free  STFD  is  obtained  by  substituting  (3)  into  (4), 

Dxx(t./)  =  A(0)DK(t,/)Aw(©)  (5) 

where  D^O,/)  is  the  TFD  matrix  of  s (t)  which  consists  of  auto-  and  cross-source  TFDs  as  its 
elements.  With  the  presence  of  noise,  which  is  uncorrelated  with  the  signals,  the  expected  values  of 
the  above  equations  yields 

£[Dx„(t,/)]  =  A(0)E[Dss(f,/)]Aw(0)  +  crl .  (6) 

In  the  above  equation,  a  is  the  noise  power,  I  is  the  identity  matrix,  and  E[.]  denotes  the  statistical 
expectation  operator. 

Equations  (5)  and  (6)  are  similar  to  the  mathematical  formula  which  has  been  commonly  used  in 
narrowband  array  processing  problems,  relating  the  source  correlation  matrix  to  the  sensor  spatial 
correlation  matrix.  However,  the  correlation  matrices  are  replaced  by  the  source  and  sensor  TFD 
matrices.  The  two  subspaces  spanned  by  the  principle  eigenvectors  of  DIX(t,/)  and  the  columns  of 
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A(0)  are,  therefore,  identical.  In  [13-15]  it  is  further  shown  that,  by  only  selecting  the  t-f  points  with 
highly  localized  signal  energy,  the  eigenvalues  and  eigenvectors  estimated  from  D „(*,/)  are  more 
robust  to  noise  than  their  counterparts  obtained  from  the  corresponding  data  covariance  matrix 
Rxx  =  £[x(t)xff(t)].  This  implies  that  key  problems  in  various  array  processing  applications  can  be 
addressed  and  solved  using  a  new  formulation  that  is  more  tuned  to  nonstationary  signal 
environments. 


3  Spatial  Polarimetric  Time-Frequency  Distributions 

3.1  Polarimetric  Time-Frequency  Distributions 

In  passive  radar,  sonar,  and  most  communication  problems,  the  received  signal  with  dual 
polarizations  can  be  expressed  as 

x(t)  =  [x[p\t)  xui(f)]r  (7) 

where  the  superscripts  (.)w  and  (.)I,!,  respectively,  denote  two  orthogonal  polarizations.  They  can  be, 
for  example,  vertical  and  horizontal  polarizations,  or  right-hand  and  left-hand  circular  polarizations. 

In  active  radar  and  sonar  applications,  the  received  signal  with  dual  transmit  and  dual  receive 
polarizations  can  be  expressed  as 

x(r)  =  [x,w,ri)  Jtl,pl(0  xlw,(r)]r  (8) 

where  the  first  letter  of  the  superscript  denotes  the  transmit  polarization,  and  the  second  letter  denotes 
the  receive  polarization.  For  notation  simplicity  and  uniformity,  we  focus  only  on  the  pp  and  qq 
components,  and  let  xlp\t)  and  xU](t)  denote  xlpp'(t)  and  x'pp,(/) ,  respectively.  In  this  way,  Eq.  (7) 
can  be  used  to  represent  both  passive  and  active  signal  processing. 

The  self-  and  cross-polarized  TFD  are  expressed  as 


141 


0./y.i(f ,/)=  -u,T)xii](t  +  r/2)(x|,|(r - r/2))’  e~n*fTdudT 


(9) 


and 

D,,y.i(f./)=  JJW* - U, t)xu] (t  +  t/2)(x[t](t - r/2))*  e~J**ftdudr  (10) 

where  the  superscripts  i  and  fc  denote  either  p  or  q.  The  self-  and  cross-polarized  TFDs  together  define 
the  polarization  TFD  matrix, 

D„(f,/)=  -  u,  T)x(t  +  T/2)x“(t  -  r/2)  e'^^dudr  .  (11) 

The  polarization  TFD  matrix  is  of  dimension  2x2,  although  it  can  be  4x4  if  the  full  four  element 
representation  in  Eq.  (8)  is  used.  The  diagonal  elements  of  D^r,/)  are  the  self-polarized  terms 

DjHjnfaf)  9  whereas  the  off-diagonal  elements  are  the  cross-polarized  terms  D^k](t,f) ,  i*k  . 

Accordingly,  the  polarization  TFD  matrix  can  be  used  to  estimate  the  self-  and  cross-polarization 
signatures  of  propagation  channels. 

3.2  Spatial  Polarimetric  Time-Frequency  Distributions 

Equations  (7)-(ll)  correspond  to  a  single  dual-polarization  sensor  case.  When  an  array  consisting 
of  m  dual-polarized  sensors  is  considered,  the  data  vector,  for  each  polarization  i,  is  expressed  as, 

xl,1(/)  =  AIO(0)s,i,(/)  +  nt'1(/) .  (12) 

It  is  noted  that,  for  structured  mixing  matrices,  they  are  polarization-independent,  i.e., 
AI,,|(©)  =  A[?1(0)  =  A(0) ,  and  Eq.  (12)  simplifies  to  the  following, 

x(,,(<)  =  A(0)sw,(f)  +  n[,1(/) .  (I3) 

The  generalization  of  Eq.  (9)  to  the  multi-sensor  receiver  is  obtained  using  Eq.  (13).  The  STFD 
matrix  introduced  in  Eq.  (4)  can  be  defined  for  each  polarization 
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D - u, r)xm (t  +  r/2)(xm (t - t/2)Y  e-Jlxftdudr . 


(14) 


In  the  noise-free  environment 

D^„(t,/)  =  A(0)Dtl,y,l(f,/)AH(0).  (15) 


In  a  similar  manner,  the  cross-polarization  STFD  matrix  between  the  data  vectors  with  two  different 
polarizations  can  be  obtained  as, 

D  ^..(t,/)  =  A(0)Ds,,y1,  (0) .  (16) 

We  are  now  in  a  position  to  tie  the  polarization,  the  spatial,  and  the  t-f  properties  of  the  signals 
incident  on  the  antenna  array.  Based  on  Eq.  (12),  the  following  vector  can  be  constructed  for  both 
polarizations. 


=B(0)s(O+n(r) 


0 

A(0)j[sl?1(oj  kj(oj 


(17) 


where  ®(®)  — 


A(0)  0  1 

0  A(0)J 


is  block-diagonal  and  s w(f),  i  =  p,q,  are  the  source  signal  vectors  for 


polarization  i.  The  STFD  of  the  dual-polarization  vector,  x(f),  can  therefore  be  defined  as 


D„(f,  /)  =  JjVri  r)x(r  +  t/2)\h  (t-r/2)  e~J2xfIdudr 

A(0)  0  I  Dsi,y,i  (f,  /)  Ds„y„  (t,  _/)!  A(0)  0  1  (18) 

=  _  0  A(0)J  (f,/)  Dsi»js0]  /)jL  »  A(0)J 


D„(/,/)  is  referred  to  as  the  spatial  polarimetric  time-frequency  distribution  (SPTFD)  matrix.  This 
distribution  serves  as  a  general  framework  within  which  typical  problems  in  array  processing  can  be 
addressed  [32-35].  In  the  next  section,  the  polarimetric  time-frequency  ESPRIT  (PTF-ESPRIT)  [32]  is 
shown  as  an  application  example  of  the  SPTFD. 
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4  Polarimetric  Time-Frequency  ESPRIT 

In  order  to  achieve  the  rotational  invariance  in  the  array  at  hand,  we  consider  a  uniform  linear 
array  with  m  cross-dipoles.  This  array  is  divided  into  two  overlapping  subarrays  of  m- 1  elements.  Let 
the  first  subarray  be  composed  of  the  left-most  m- 1  cross-polarized  antennas,  and  the  second  subarray 
be  composed  of  the  right-most  m- 1  cross-dipoles.  Additionally,  let  the  array  response  matrices  of  the 
identically  polarized  sensors  of  the  two  subarrays  be  A,(0)  and  A  2  (0) ,  respectively.  Accordingly, 

Aa<0)  0  1  jA.f©)  0  TO.  01 

o  a2(0)'j  [  o  A,(@yj[o  ; 


where  the  rotation  operator  <&  can  be  modeled  as 


<t>  =  diag[e 


in(#, )  -;2JTySin(«, )  -j2xjsin(8, ) 

A  ,e  A  A  ], 


d  is  the  interelement  spacing,  and  X  is  the  wavelength  of  the  impinging  signals. 

By  performing  joint  block-diagonalization  [12,36]  on  the  SPTFD  matrices  Dxx(/,/>  over  multiple 
(r,/)  points  where  the  energy  of  the  signal  arrivals  is  concentrated,  we  obtain  the  signal  and  noise 
subspaces,  represented  as  matrices  Us  and  U„ ,  respectively  . 

The  signal  eigenvectors,  comprising  the  columns  of  U, ,  span  approximately  the  signal  subspace  of 
dimension  n  such  that  there  exists  a  transformation  matrix  T  that  satisfies 


0  |]T 

A(0)J 


(21) 


where  A(0)  is  the  full  (undivided)  array  response  matrix  for  one  polarization. 

By  applying  the  same  transformation  matrix  T  to  the  steering  matrices  of  the  two  subarrays,  the 
following  matrices  are  defined, 
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polarization  parameters,  yand  <p.  representing  the  amplitude  ratio  (tan0— and  phase 

difference  between  the  two  polarization  components,  respectively.  The  normalized  frequency  of  the 
sinusoid  is  0.1.  All  signals  have  the  same  signal  power  (SNR=5  dB).  The  data  length  is  256  samples. 

The  pseudo  Wigner-Ville  distribution  (PWVD)  is  used  to  compute  the  TFDs  where  a  rectangular 
window  of  length  65  is  used. 


Table  2.  Signal  parameters 


DOA 

(deg) 

Y (deg) 

(deg) 

Source  1 

-3 

45 

0 

Source  2 

3 

45 

180 

Source  3 

5 

20 

0 

5.1  Array  and  Polarization  Averaging  of  TFDs 

The  proper  selection  of  t-f  autoterms  is  often  important  in  array  processing  based  on  STFDs  and 
SPTFDs  [14,16,38].  The  presence  of  crossterms  and  noise  often  obscure  the  identification  and. 
selection  of  t-f  autoterm  regions.  Averaging  TFDs  across  different  sensors  and  different  polarizations 
can  suppress  the  effect  of  noise  and  crossterms  and  render  it  easier  to  identify  the  autoterm  TFDs 
[39,40].  The  suppression  of  crossterms  is  highly  dependent  on  the  spatial  correlations  and  polarization 
diversity  among  the  signals. 

Figure  1  shows  the  PWVD  of  the  signal  received  at  the  first  vertical  sensor.  With  the  single  sensor 
and  single  polarization,  the  crossterms  between  the  source  signals  and  the  noise  obscure  the  correct 
identification  of  the  autoterm  of  each  signal  component.  Figure  2  shows  the  PWVD  averaged  across 
the  four  vertical  sensors.  While  the  noise  is  substantially  suppressed,  the  array  averaging  does  not 
significantly  reduce  the  crossterms  because  of  the  close  orientation  of  the  sources.  To  further  suppress 
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the  crossterms,  we  utilize  both  the  spatial  and  polarization  information.  Figure  3  shows  the  PWVD 
averaged  over  the  four  sensors  and  both  polarizations.  In  this  example,  because  source  1  and  source  2 
have  orthogonal  polarizations,  their  associated  crossterms  are  completely  suppressed. 


time 


Figure  1.  PWVD  of  the  vertical  component  received  at  the  first  antenna  sensor. 


147 


0| 


0.1  v  ■t££A# 


fr0.2[ 

§ 

Jo.3 

0.4- 


,.A  ,.:  • 
■  rv  <r< 


0.5!- 


50  100  150  200  250 

time 


Figure  3.  PWVD  averaged  across  four  antennas  and  both  polarizations. 


5.2  PTF-ESPRIT 

For  both  time-frequency  ESPRIT  (t-f  ESPRIT)  [18]  and  PTF-ESPRIT,  the  array  averaged  PWVD 
is  first  used  to  identify  the  autoterm  regions.  The  search-free  PTF-ESPRIT  provides  a  DOA 
estimation  which  is  compared  to  that  of  the  conventional  ESPRIT,  polarimetric  ESPRIT,  and  t-f 
ESPRIT.  In  the  underlying  situation  where  the  source  signatures  can  be  separated  in  the  t-f  domain, 
only  the  t-f  points  on  the  signal  signature  of  a  single  source  are  considered  for  STFD  and  SPTFD 
matrix  construction.  The  PTF-ESPRIT  outperforms  other  ESPRIT-based  methods  by  taking 
advantages  of  such  source  discriminatory  capability,  in  addition  to  the  SNR  enhancement  and 
polarimetric  selection. 

Figure  4  shows  the  root  mean  square  error  (RMSE)  performance  of  the  PTF-ESPRIT  and  other 
ESPRIT-based  methods  versus  the  input  SNR,  where  the  least-squares  approach  is  used  for  all 
methods  and  the  results  are  obtained  from  100  independent  trials.  For  conventional  and  t-f  ESPRIT, 
only  the  vertical  polarization  components  of  the  source  signals  are  used.  For  the  t-f  and  PTF-ESPRIT 
methods,  only  the  first  source  signal  is  selected  in  STFD  and  SPTFD  matrix  construction, 
respectively.  The  RMSE  performance  is  evaluated  for  the  first  source  signal.  It  is  evident  that  the 
PTF-ESPRIT  outperforms  all  the  other  methods.  It  is  clear  that  the  polarimetric  ESPRIT  provides 
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satisfactory  DOA  estimation  only  when  the  input  SNR  is  moderate,  and  the  conventional  ESPRIT 
fails  to  do  so  for  all  the  input  SNR  values  simulated.  To  the  contrary,  the  PTF-ESPRIT  provides  1- 
degree  RMSE  when  the  input  SNR  is  at  a  low  level  of  about  -7  dB. 


Figure  4.  RMSE  performance  of  ESPRIT  algorithms  versus  input  SNR. 

In  Figure  5,  DOA  estimates  of  30  independent  trials  are  shown  for  the  different  ESPRIT  methods 
utilizing  the  least-squares  approach,  where  the  input  SNR  is  5  dB.  It  is  evident  that,  at  this  low  input 
SNR  level,  only  the  PTF-ESPRIT  produces  accurate  and  consistent  estimates  of  the  DOAs  of  all  three 
signals.  While  the  t-f  ESPRIT  provides  comparable  DOA  estimation  for  source  signal  1,  the  variance 
of  the  DOA  estimations  is  much  greater  for  the  other  two  signals. 
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Conventional  ESPRIT  Polarimetric  ESPRIT 


TF  ESPRIT  PTF-ESPRIT 


Figure  5.  *  DOA  estimation  results  from  different  ESPRIT  methods. 

6  Conclusion 

The  concept  of  SPTFD  has  been  proposed  and  shown  to  be  a  powerful  platform  to  utilize  the 
polarization  and  temporal  signatures  of  signal  arrivals  for  sophisticated  array  processing.  To 
demonstrate  the  advantage  of  the  SPTFD  platform,  we  have  considered  the  DOA  estimation  problem 
and  proposed  the  polarimetric  time-frequency  ESPRIT  (PTF-ESPRIT)  method  as  an  example  of  its 
applications.  In  such  scenarios  where  the  signals  are  highly  localized  in  the  time-frequency  domain 
and  diversely  polarized,  the  proposed  PTF-ESPRIT  significantly  outperforms  the  other  existing 


150 


ESPRIT  methods,  including  conventional  ESPRIT,  time-frequency  ESPRIT,  and  polarimetric 
ESPRIT. 
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